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ABSTRACT 


The  dynamic  stall  of  rapidly  pitching  and  oscillating 
airfoils  is  investigated  by  the  numerical  solution  of  the 
full  compressible  unsteady  l /.o-dimensionai  Navier-btoKes 
equations  using  an  alternating-direction-implicit  scheme. 

The  flow  is  assumed  to  be  fully  turbulent,  and  the  turbulent 
stresses  are  modelled  by  the  Baldwin-Lomax  eddy  viscosity 
mode].  Three  airfoils  (NACA  0012,  NACA  0012-33,  and  NACA 
0012-63)  are  analyzed  for  the  purpose  of  examining  the 
influence  of  leading-edge  geometry  on  unsteady  flow  separa¬ 
tion.  It  is  found  that  a  larger  leading  edge  radius,  thick¬ 
er  contouring  of  the  forward  part  of  the  airfoil,  or  in¬ 
creasing  reduced  frequency  results  in  delaying  flow  separa¬ 
tion  and  formation  of  the  dynamic  stall  vortex  to  a  higher 
angle  of  attack,  yielding  higher  peak  C-.  Within  the  scope 
of  this  study,  the  pressure  gradient  encountered  by  the  flow 
at  initial  separation  is  found  to  be  independent  of  reduced 
frequency  and  freestream  speed.  The  critical  pressure 
gradient  is  dependent  on  leading  edge  radius  and  increases 
for  decreasing  leading  edge  radius. 
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INTRODUCTION 


I . 

A.  BACKGROUND 

The  aeronautical  community  has  been  aware  of  the  phenome¬ 
non  of  dynamic  stall  for  several  decades.  Dynamic  stall  is 
characterized  b-,  separated  flow  and  shedding  of  the  leading 
edge  vortex  from  the  upper  surface  of  an  airfoil  which  is 
rapidly  pitched  to  angles  of  attack  beyond  the  normal  stalling 
angle  of  attack.  The  phenomenon  results  in  significant 
temporary  increases  and  decreases  in  lift,  drag,  and  moment 
coefficients . 

The  scope  of  this  work  is  the  investigation  of  the 
phenomenon  of  dynamic  stall  for  rapidly  pitching  and  oscillat¬ 
ing  airfoils.  The  phenomenon  of  dynamic  stall  was  observed 
for  the  first  time  for  flows  over  the  retreating  blades  of  a 
helicopter.  The  dynamic  stall  resulting  from  the  oscillatory 
motion  of  the  rotor  blade  is  associated  with  an  increase  in 
lift  and  the  development  of  a  severe  nose  down  pitching 
moment.  The  effects  of  dynamic  stall  are  usually  undesirable 
for  the  helicopter,  and  where  possible  special  care  is  taken 
to  reduce  its  effects  by  special  design  of  the  rotor.  On  the 
other  hand,  interest  has  recently  developed  to  exploit  the 
increased  lift  obtained  from  the  rapid  pitch-up  motion  of  an 
airfoil  in  order  to  enhance  the  maneuverability  and  extend  the 
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flight  envelope  of  the  modern  fighter  aircraft  to  the  high 
angle  of  attack  regime,  or  to  alleviate  retreating  blade  stall 
for  helicopters  utilizing  Higher  Harmonic  Control  (HHC). 

It  is  well  known  that  for  flow  over  an  airfoil  at  fixed 
angle  of  attack  the  streamlined  airflow  is  disrupted  once  a 
critical  angle  of  attack  is  exceeded.  At  stall,  the  flow  over 
the  upper  surface  of  the  airfoil  separates  and  the  lift  drops. 
It  was  observed,  however,  that  rapid  pitch-up  motion  of  the 
airfoil  delays  static  stall  so  that  high  lift  tan  be  main¬ 
tained  for  angles  of  attack  beyond  the  static  stall  angle. 

As  the  pitch-up  angle  exceeds  the  static  stalling  angle  of 
attack,  a  thin  layer  of  reversed  flow  develops  in  the  boundary 
layer.  This  reversed  flow  occurs  in  two  types  of  stall: 
trailing  edge  stall  and  leading  edge  stall  [Ref.  1J.  In 
trailing  edge  stall,  the  reverse  flow  region  begins  near  the 
trailing  edge  and  traverses  forward;  in  leading  edge  stall, 
the  reverse  flow  region  occurs  near  the  suction  peak  just  aft 
of  the  leading  edge.  In  both  cases,  a  vortex  begins  to  form 
near  the  leading  edge  region,  expands,  and  moves  downstream. 
The  angle  of  attack  at  which  the  vortex  is  formed  depends  on 
airfoil  shape,  pitch  rate,  mean  angle  and  amplitude,  Mach 
number,  and  Reynolds  number. 

The  vortex  formed  at  the  leading  edge  is  called  the 
dynamic  stall  vortex,  and  moves  with  a  speed  of  approximately 
0.4  freestream  speed  relative  to  the  airfoil  as  the  pitch-up 
progresses.  Lift,  drag,  and  pitching  moment  increase  signifi- 
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cantly  until  the  vortex  approaches  the  trailing  edge,  then 
drop  sharply,  but  not  simultaneously  (Figure  1).  The  unsteady 
surface  pressure  increases,  and  the  suction  peak  appears  at 
different  locations  along  the  chord  as  the  dynamic  stall 
vortex  moves  over  the  airfoil.  Secondary  and  tertiary 
vortices  may  also  be  present  and  produce  additional  suction 
peaks  and  fluctuations  in  airloads. 

To  date,  most  theo¬ 
retical  investigations 
on  dynamic  stall  have 
focused  on  the  effects 
of  variations  of  reduced 
frequency  (k  =  &c/2UJ  , 
angle  of  attack,  free 
stream  speed,  and 
Reynolds  number  on  a 
particular  airfoil  shape 
(usually  the  NACA  0012 
airfoil).  McCroskey 

[Ref.  1]  documented  the 
effects  of  reduced  fre¬ 
quency,  amplitude  and 
Mach  number  on  different 
airfoils.  These  experi¬ 
mental  studies  demonstrated  the  large  hysteresis  in  the 


Figure  1.  Events  of  Dynamic  Stall 
on  NACA  0012  Airfoil 
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occurrence  of  stall  of  an  oscillating  airfoil  compared  to  a 
fixed  angle  of  attack  airfoil.  McCroskey  and  Pucci  [Ref.  2] 
identified  varying  regimes  of  viscous-inviscid  interaction 
during  varying  degrees  of  unsteady  flow  separation.  The 
conclusion  of  the  experimental  studies  of  Refs.  1  and  2  was 
that  the  reduced  frequency  has  a  dominant  effect  on  the 
development  and  progression  of  the  dynamic  stall.  Experi¬ 
mental  work  by  Chandrasekhara  and  Carr  [Ref.  3]  using  the  NACA 
0012  airfoil  showed  that  a  dynamic  stall  vortex  always  forms 
near  the  leading  edge  of  an  oscillating  airfoil.  Their  study 
also  documents  the  movement  of  the  dynamic  stall  vortex. 
Chandrasekhara  and  Brydges  [Ref.  4]  documented  the  effects  of 
increasing  amplitude  on  an  oscillating  airfoil  in  both 
compressible  and  incompressible  flow.  They  showed  that  larger 
amplitudes  resulted  in  vortex  retention  at  higher  angles  of 
attack  for  a  given  Mach  number  and  reduced  frequency. 

Progress  in  computational  fluid  dynamics  has  made  possible 
the  study  of  dynamic  stall  by  numerical  solution  of  the 
unsteady  Navier-Stokes  equations.  Mehta  [Ref.  5]  demonstrated 
that  Navier-Stokes  simulation  of  the  unsteady  incompressible 
flow  around  the  airfoil  in  oscillatory  motion  can  reproduce 
the  experimentally  observed  results.  Wu  et  al.  [Ref.  6] 
presented  solution  procedures  based  on  an  integral  formulation 
of  the  incompressible  Navier-Stokes  equations  for  the 
computation  of  unsteady  flow  over  airfoils.  Beddoes  [Ref.  7] 
and  Jang  et  al.  [Ref.  8]  presented  viscous-inviscid  computa- 
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tion  methods  for  unsteady  flows.  Sankar  and  Tang  [Ref.  9], 
Visbal  [Ref.  10],  and  Ekater inaris  [Ref.  11]  used  ADI  (Alter¬ 
nating  Direction  Implicit)  numerical  schemes  for  the  solution 
of  the  compressible  Navier-Stokes  equations,  and  investigated 
the  effects  of  compref  ibility  on  dynamic  stall.  The  numeri¬ 
cal  solutions  for  both  incompressible  and  compressible  flows 
showed  good  agreement  with  experimental  results,  and  predicted 
the  events  of  dynamic  stalls. 

B.  PURPOSE 

As  indicated,  the  aforementioned  studies  focused  on  the 
variation  of  flowfield  parameters  on  the  dynamic  stall  of  a 
particular  airfoil.  The  purpose  of  this  study  is  the  system¬ 
atic  investigation  of  the  effect  of  leading  edge  geometry  on 
the  development  of  dynamic  stall.  It  is  expected  that  the 
effect  of  the  leading  edge  geometry  will  have  a  significant 
effect  on  both  the  development  of  the  dynamic  stall  vortex  and 
its  subsequent  shedding.  This  numerical  investigation  is 
intended  to  provide  a  cost-effective  means  of  quantifying 
which  airfoil  parameter  variations  should  be  analyzed  in  more 
expensive  wind-tunnel  tests.  The  numerical  investigation  of 
airfoil  parameter  variations  offers  the  benefit  of  optimizing 
the  utilization  of  more  costly  test  facilities. 

The  investigation  is  conducted  using  a  numerical  solution 
of  the  full  two-dimensional,  Reynolds-averaged  Navier-Stokes 
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Equations,  and  the  Baldwin-Lomax  eddy  viscosity  model  of  Ref. 
16  is  used  to  obtain  the  turbulent  stresses. 

C.  AIRFOIL  SELECTION 

Modif icat ions  to  the  NACA  0012  airfoil  were  made  forward 
of  the  point  of  maximum  thickness  (12%  thick  at  30%)  chord  by 
the  method  suggested  in  Ref.  12.  The  NACA  0012  profile  was 
retained  aft  of  30%  chord  for  the  two  new  airfoils.  In  this 
manner,  changes  in  the  flow  could  be  attributed  directly  to 
the  changes  in  the  forward  part  of  the  airfoil.  The  two  air¬ 
foils  consisted  of  the  NACA  0012-63  and  NACA  0012-33  airfoil 
section  forward  of  30%  chord.  The  NACA  0012-63  has  the  same 
leading  edge  radius  as  the  NACA  0012  (1.58%  chord)  but 

different  curvature  to  the  point  of  maximum  thickness.  The 
NACA  0012-33  has  a  smaller  leading  edge  radius  (0.39%  chord) 
with  curvature  necessary  to  achieve  the  same  12%  thickness 
ratio  at  30%  chord.  Airfoil  coordinates  are  provided  in  Table 
1,  and  the  resulting  modified  profiles  are  shown  in  Figure  2. 
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TABLE  1 


AIRFOIL  COORDINATES 


x/c 

N0012 

.01 

.01701 

.  02 

.02360 

.03 

.  02842 

.  04 

. 03232 

.05 

.  03555 

.  06 

. 03840 

.  07 

.04090 

.08 

.  04313 

.  09 

. 04505 

.  10 

. 04681 

.11 

. 04842 

.12 

. 04990 

.  13 

. 05123 

.  14 

. 05242 

.15 

.  05345 

.  16 

. 05441 

.  18 

.05610 

.  20 

.05742 

.22 

.  05840 

.  24 

.05903 

.  25 

. 05943 

.  26 

. 05964 

.  28 

.  05992 

.  30 

. 06000 

.  40 

.  05803 

.  50 

. 05294 

.  60 

.  04563 

.  70 

. 03664 

.  80 

.02623 

.  90 

.01448 

.  95 

.  00807 

1 . 00 

.  00126 

N0012-63  N001 2-33 

.01682  .01040 
.02320  .01550 
.02783  .01969 
.03160  .02332 
.03474  .02661 
.03750  .02962 
.03991  .03241 
.04210  .03500 
.04403  .03741 
.04580  .03966 
.04742  .04177 
.04890  .04374 
.05023  .04557 
.05147  .04734 
.05260  .04888 
.05362  .05035 
.05540  .05300 
.05686  .05515 
.05802  .05692 
.05890  .05830 
.05924  .05801 
.05952  .05925 
.05988  .05982 
.06000  .06000 
.05803  .05803 
.05294  .05294 
.04563  .04563 
.03664  .03664 
.02623  .02623 
.01448  .01448 
.00807  .00807 
.00126  .00126 


Figure  2.  Airfoil  Profile  Comparison  • 
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II.  GOVERNING  EQUATIONS 


The  flow  of  a  compressible,  viscous  fluid  satisfies 
conservation  of  mass,  momentum,  and  energy.  The  conservation 
of  mass  is  expressed  by  the  Continuity  Equation,  the  conserva¬ 
tion  of  momentum  by  the  Navier-Stokes  Equations,  and  the 
conservation  of  energy  by  the  Energy  Equation.  The  conserva¬ 
tion  equations  are  derived  in  the  following  section.  The 
derivation  is  done  in  an  Eulerian  frame  of  reference. 

A.  CONTINUITY  EQUATION 

Consider  a  control  volume  in  a  flow  field  where  flow 
properties  vary  with  both  time  and  space.  Conservation  of 
mass  requires  that  the  rate  of  change  of  mass  inside  the 
control  volume  equals  net  mass  flux  out  of  the  control  volume. 
This  is  expressed  as: 

-£///„<,,<’  d(voJ)  *  /,/p|7d3  - 0  11  > 

using  Gauss's  theorem  for  a  surface  integral,  Eq .  1  becomes 

//Llfd(TOi)  +  ffLv  (r*  d(vol>  ■  0  (2> 

or , 

!!Llt  * ]  d ( vol )  =  0  ( 3 ) 

therefore,  for  an  arbitrary  control  volume, 


9 


0 


(4) 


|f  ^  V-(p^  - 

for  any  continuous  flow.  For  a  two-dimensional  flow,  and 
Cartesian  coordinates,  Eq .  4  reduces  to 


ftp  +  ft(pu)  +  ft(py)  _  0 

dt  dx  dy 


(5) 


B.  THE  MOMENTUM  { NAV I ER- STOKES )  EQUATIONS 

Consider  a.  fluid  element  in  a  rectangular  Cartesian 
coordinate  system.  The  stresses  and  pressures  are  shown  in 
Figure  3,  and  body  forces  are  neglected.  Summation  of  forces 
in  the  x-direction  yields: 


Fx  =  pdxdydz-^ 

=  (p+ J  dydz  +  f  (t  ^  •»  dx)  -t  ^3  dydz 

*  [<,xyx+^dy)-x„\dxdz  *  [t  tx*^fdz)-xzx\dxdy 


(6) 


or  , 


Du  _  _  _ftp  + 


ftt 


dx 


yx 


dt 


(7) 


Dt  dx  dx  dy  dz 
Similarly,  summation  of  forces  in  the  y-  and  z-  directions 
yields: 


Dv  _  _  dp_  + 


OT 


xy 


OX 


Dt 


dy 


dx 


yy 


ox 


dy 


ty 


dz 


(8) 


Dw  =  _dp  + 


Dt 


dz 


d'  xz  +  y*  +  dx  rr 

dx  dy  dz 


(9) 
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Figure  3.  Forces  on  a  Fluid  Element 


Using  the  Continuity  Equation,  the  equations  expressing  the 
conservation  of  momentum  for  a  two-dimensional  flow  are: 


^  *4  «”'■**>  -  ax 


dx 


xy 


dy 


(10) 


+  a  (Puv)  +  /(Pvj  +  P)  =  -2^  +  ilsz  (id 

dt  dx  oy  ox  dy 

The  relation  of  the  viscous  stresses  tvv,  and  r  ,  to  the 

**  yy  xy 

independent  variables  is  developed  in  the  following  discus¬ 
sion. 

For  a  two-dimensional  flow  in  Cartesian  coordinates  with 
an  infinitesimal  fluid  element  undergoing  distortion  due  to 
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stresses  as  shown  in  Figure  4,  the  angular  displacements  A0, 
and  A0;  are: 

49,.  -git  49,.-g4e  <12> 

The  strain  increment  is  given  by: 

A K  =  A02  -  A0,  ,  <13) 

in  the  limit,  the  rate  of  strain  is  given  by: 

dk  =  ^Z  +  iH  (14) 

1  xy  dt  dt  dt  dx  dy 

Using  Newton’s  Law  of  Fluid  Friction  (definition  of 
viscosity) , 


xt)  = 


(15) 
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then  the  rate  of  strain  caused  by  the  tangential  shear  stress 
is  given  by: 


For  large  velocity  gradients  the  normal  stresses  t.  and  t  . 
caii  be  significant  and  result  in  a  viscous- induced  normal 
force  on  the  fluid  element.  For  example,  as  documented  by 
Schlichting  [Ref.  13],  in  order  for  fluid  isotropy  to  be 
maintained  at  every  point,  the  principal  axes  of  stress  and 
rate-of-strain  must  coincide  to  avoid  introducing  a  preferred 
rotation  direction.  With  this  concept,  a  normal  stress  must 
depend  both  on  its  respective  component  rate  of  strain  as  well 
as  the  shearing  strain  rates,  with  different  weighting 
factors.  Choosing  the  factor  2p  for  the  component  direction 
factor,  which  causes  Newton's  Law  of  Friction  to  be  satisfied 
for  simple  shear,  obtain 


,  4b)  *  2H|h 

ox  oy  ex 


where  X  is  the  shearing  stress  proportionality  factor.  Using 
Stokes's  hypothesis,  X  =  -2/3  p,  obtain 


(19) 


C.  THE  ENERGY  EQUATION 

Conservation  of  energy  is  the  manifestation  of  the  First 
Law  of  Thermodynamics  (dE  =  dW  +  dQ )  to  a  moving  fluid  element 
in  rectangular  Cartesian  coordinates.  The  First  I.aw  of 
Thermodynamics  is  applied  to  a  control  volume  (Figure  5), 
where  the  energy  fluxes  are  shown.  The  rate  of  change  of 
energy  inside  the  fluid  element  is  equal  to  the  net  flux  of 
heat  into  the  element  plus  the  rate  of  work  done  on  the 
element  by  pressure  and  viscous  stress. 

The  rate  of  change  of  energy  of  the  fluid  element  having 
an  instantaneous  internal  energy  per  unit  mass  e  and  speed  V 
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[gx  -  (<jx  +  dx)  ]  dydz  -  -  dxdydz  (  21 ) 

Accounting  for  the  y  and  z  directions,  the  total  heat  trans¬ 
ferred  out  of  the  fluid  element  is: 


+  ^~)  dxdydz  =  -  (V-<^)  dxdydz  (22) 

dx  dy  dz 

where  the  heat  flux  is  expressed  in  terms  of  the  temperature 
gradient  according  to  Fourier's  law  of  heat  conduction  as 


<ki  = 


dT 

dx± 


(23) 


and  k  is  the  thermal  conductivity,  considered  constant  and 
independent  of  temperature. 

The  rate  of  work  done  on  the  fluid  element  in  the 
x-direction  by  pressure  and  shear  stresses  is 
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r  3  (UP)  +  ^(UTxx)  +  ^UTyx> 
dx  dx  dy 


^  —  ]  dxdydz  (24) 

dz 


+ 

Accounting  for  the  y  and  z  direction  pressure  and  shear 
stresses,  the  net  rate  of  work  done  on  the  fluid  element  can 
be  expressed  as: 


.  _  -  8(ut„J  d(uxvx )  d  ( ux  ) 

[-V*pV  + - _  -iE£  + - _>2L_  + 


8x 


0y 


8z 


divXxy)  3(vtvv)  d{vxzv) 


yy, 


dx 


dy 


zy 


dz 


(25) 


d(wx  xz)  d(m  yz)  d{wx  zz) 


}  dx  dy  dz 


dx  dy  dz 

The  complete  energy  equation  for  a  viscous  flow  with  no 


external  heating  can  then  be  expressed  as: 


p-£r(e+4^2)  =  t-V*g-  V-pV 

(it  2 


.  diux^)  _  d(uxyx)  _  3(ut,x) 
8x  3y  3z 

,  ,  SiVXyy)  t  d(vx  zy) 

dx  dy  dz 


,  d(^xz)  ,  d(^yZ> 


3x 


8y 


3z 


]  dxdydz 


The  two-dimensional  compressible  viscous  flow  energy  equation 
reduces  to : 


P—  (e+±V2)  4 -A  (pu)  +^-(pv) 
dt  2  dx  dy 


Qx  ^  ^  xx  +  ^  xy  ^x^  +  0y  (  yx+  ^  w 


(27) 


yy 


Using  the  continuity  equation  and  setting 
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the  total  energy  per  unit  volume,  the  energy  equation  becomes 


BE 

dt 


+ 


~  [  (F+p)  u]*4-[  {E+p) 

Bx  ay 


v ) 


4-iu< 

dx 


+  VT 

xx  xy 


d_ 

dy 


4x>  +^-(ut 


+  VT 

xy  yy 


-  <?v) 


(29) 


where  t  t  and  t  vu  are  as  previously  derived  (see  Eq .  19). 

AA  >y  Ay 

D.  CONSERVATION  LAW  FORM  OF  THE  GOVERNING  EQUATIONS 

The  conservation  law  form  of  the  governing  equations  for 
the  two-dimensional  viscous  compressible  flows  can  be  written 
as  [Ref.  14] 


dQ  ^  dF  x  dG  BR  dS  ,q0, 

if  =  1  ’ 

where  Q  is  the  vector  of  dependent  variables,  F  and  G  are  the 
inviscid  fluxes  along  the  x  and  y  directions,  respectively, 
and  R  and  S  are  the  viscous  fluxes  along  the  x  and  y  direc¬ 
tions,  respectively.  These  terms  are  given  by: 


~p ' 

PU 

p  V 

pu 

F  = 

pu2+p 

G  = 

puv 

pv 

puv 

p  v2+p 

_E 

jE+p)  u 

( E+p)  v 

—  — 

—  -H 

0 

0 

T  XX 

S  = 

^  xy 

^  xy 

t  yy 

xx+V^  X3  ~<?x 

UT.'0'+VTyy-<7y 

Non-dimensionalizing  the  equations  with 


P*  = 


P 


T’  = 


T_ 

Tm 


f  = 


t 

L/Vm 


where  L  is  the  reference  length,  and  setting 


Re 


P  mVmL 

P. 


(31) 


the  Reynolds  number  based  on  the  reference  length,  obtain  the 
non-dimensional  form  of  the  conservation  law: 

+  +  =  _±_  ilE  +  JH)  02) 

dt  dx  dy  Re  dx  dy 

with  non-dimensional  variable  Q  and  the  inviscid  flux  terms  (F 
and  G)  as  before,  and  the  viscous  terms  given  by: 


R  = 


xy 


U  T  ■+  VT  - 
xx  xy 


n  dr 

(Y-l)  M2Pr  dx 


(33) 
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0 


s  = 


(34) 


L  (y-1)  M2Pr  dy\ 

Here  Pr  -  Cp  p/K  is  the  Prandtl  number,  and  M  =  jV|/aw.  In 
order  to  enable  solution  of  the  governing  equation  for 
arbitrary  geometries,  the  governing  equations  are  expressed 
for  a  curvilinear  generalized  coordinate  system.  The 
curvilinear  coordinate  system  (  (*  ,  T) )  is  linked  to  the  Cartesian 
coordinate  system  (x,y)  by 


I \  =  \{x,y,t)  ,  t|=ti(x, y,t)  (35) 

The  curvilinear  coordinate  space  (  £  ,  ^  )  is  referred  to  as  the 
computational  domain  and  is  linked  to  the  physical  domain 
(x,y)  by  the  non-zero  Jacobian  of  the  coordinates  transforma¬ 
tion 


j  =  iiidii  =  ilia,  ilia  = _ i _ 

3{x,y )  dx  dy  dy  dx  dx  j)y  _  Bx^  By  (36) 

81;  3t)  8t]  5^ 

It  can  be  shown  that  the  governing  equations  [Ref.  18]  retain 
the  strong  conservation  law  form  for  a  generalized  coordinate 
system.  The  two-dimensional  strong  conservation  law  form  for 
a  generalized  coordinate  system  is: 

dQ  .  dF  dG  _  1  /  dR  BS  *  { 37 ) 

Bt  31  at]  "  Re  v  as  dr]  1 
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where 
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0 


[  (ut 
+ 


t  il+T 
**  dx  *  dy 


in 


T 


in 

xy  dx 


+  T 


in 

^  dy 


+  VT 

xy  yy 


(UT  ^y+VT 


_ ^ _ dT  j  _drj_ 

(Y-l )M2Pr  dy  dx 

--H _ ii)  inj 

yy  (y-l)  m2  Pr  dy  dy 


(44) 
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Ill . 


NUMERICAL  APPROACH 


A.  NUMERICAL  PROCEDURE 

The  strong  conservation  law  form  of  the  two-dimensional 
Continuity,  Navier-Stokes ,  and  Energy  Equations  in  generalized 
coordinates  presented  in  Chapter  II  provides  a  useable  format 
for  implementing  a  numerical  solution  technique.  The  numeri¬ 
cal  method  used  for  the  integration  of  the  governing  equations 
is  a  finite  difference  numerical  scheme  based  on  the  Beam- 
Warming  algorithm  [Ref.  15].  The  viscous  terms  are  retained 
in  both  directions  in  order  to  enable  capturing  intense 
viscous  effects  encountered  in  massively  separated  flow 
regions  at  high  angles  of  attack.  The  turbulent  stresses  were 
modeled  using  the  Baldwin-Lomax  eddy  viscosity  model  [Ref. 
16].  An  algebraic  C-type  grid  (157x58)  was  used  for  the 
computations.  The  boundary  conditions  were  treated  explicit¬ 
ly,  and  the  unsteadiness  was  imposed  by  the  motion  of  the 
grid. 

B.  THE  BEAM-WARMING  ALGORITHM 

The  strong  conservation  law  form  of  the  two-dimensional 
compressible  Navier-Stokes  equations  in  the  vector  notation  is 
as  follows: 
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(45) 
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where  An  U  =  Un+1  -  Un.  By  substituting  Eq .  (45)  into  Eq .  (47)  , 


obtain 


E  4-  A”Vj  +  A"Vj)+  (— ArlF  +  A 


4"u  ■'T+T;  [s  <-a"E  +  a"vi  +  A"V,)+  ~  (-a"F  +  a"v 

+  ttt,  [s (~c” + v" + v;> + £  <-F" +w"+ 'v;)l 

+  r+97a""'u  +  0  (A0’+  (Ai)’l 


"W,  4-  AnWj) 


(48) 


The  delta  terms  are  linearized  using  truncated  Taylor  series 
expansions,  so  that: 

En*l  =  En+  (_|l)  a  (yn^-jyn)  +0[(At)2]  (49) 

dU 

which  can  be  rewritten  as 

A  nE  =  [A]  nAnt7+0[  (At)2]  (50) 

where  [A]  is  the  flux  Jacobian  matrix  3E/dU  given  by 


\A)  =- 


3 -7  2  , i ~y  2 

+-—U1 


■+  (I  -7 )u(u7  4-  uJ) 


(T  —  3)u 


7^/  ,  7“  1  ,,  ,  , 

■— +  -y-(3“  +  y  ) 


(7  -  1  )w 


(7-  l)»i> 


0  ~7) 


In  a  like  manner,  An  F  can  be  linearized  as 


A  nF  =  [5]  nAnU+0[  (At)2] 


where  [ B]  is  the  Jacobian  matrix  dF/dU  given  by 
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[B]  =- 


0 

uv 

3 -y  7  .  1  ~y  2 
—v+—u 


yE,u 


+  (l  —  y)v(u7  +  v1) 


0 

—v 


(7-1)“ 
(y  —  1  )uv 


-1 

—u 

(7  -  3)u 

1 


yE,  7- 

P  2 


(3d1  +  ua) 


0 

0 

1  -y 

-yv 


(53) 


The  viscous  delta  term  An  V1(U,UX)  is  linearized  by  writing 


(54) 


'"v‘  =  +  (^r)" A,,u- +  0|(A,)!  1 

=  [/>]”  AnU  +  [/?]"A”UJC  +  0[(A/)2] 

=  m  -  1**1  T  A"U  +  ~  ([/?)"  A”U)  +  Ol(A02] 


where  [  P]  is  the  Jacobian  dV^dU,  [ /?]  is  the  Jacobian  dV1/dUx, 
and  [/?v  ]  =  d[R]/dx.  These  matrices  can  be  written  as 


0 

1 

1 

0  ! 

0 

!  o 

[P\- 

-  [K*]  =  - 1 

RR 

1 

I 

1 

1 

(R I 

1 

0 

1 

1  0 

1 

p 

1 

0  | 

Px 

1  0 

~u7{b)x~v1^ 

1 

1 

I 

1 

*(R  | 

Wx 

1 

!  ° 

0 

i  0 

\ 

0 

1  0 

_ 1 

4 

-j»u 

|  4 

!  3  *• 

1 

1 

1 

0 

1 

1  0 

p 

-UV 

1 

1  0 
| 

1 

1 

u 

1 

1  0 

R-£> 

-(p -A 

!k 

p 

1  /4  *' 

!V3"-r 

)« ! 

1 

(- 

k\ 
—  u 
cv/ 

1 

1  k 

J  <"u 

(55) 


(56) 
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The  matrix  for  [  P)  -  [  J?x]  is  obtained  by  assuming  that  y  and  k 
are  locally  independent  of  U.  In  a  like  manner,  An  ^(U.U^,)  is 
linearized  as 

A  nW2  =  (  [Q]  ~  [Sy]  )  nAnC7+-|;  (  [S]  nAny)  +0[  (At)2]  (57) 


where 


0  | 

0  ! 

0  | 

0 

~UHy 

Vy  | 

0  ! 

0 

lej  -  is,]  =-i 

1 

1 

0  I 

1 

| 

(R  1 

0 

-v’( 

^4  \  2  ! 

K~3  /y  ~u  11)1  i 

1 

Ufly 

1 

u(H  I 

0 

and 

0 

i 

0 

|  o 

-fJU 

i 

! 

M 

!  ° 

1 

151  “7 

- 

4 

1 

1 

1 

1 

0 

;  4 

I 

| 

iV-1- !  1 

cv)  cv  p  j 

f  *> 

\  cv/ 

)• 

0 

0 

0 

k_ 

cv 


(58) 


(59) 


The  cross-derivative  terms  can  be  evaluated  explicitly  without 
loss  of  accuracy  by  noting  that 
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A  V  =  A^Vj  +  OKAt)2] 
AnWn  =  A^^  +  Ot  (At)2] 


(60) 


for  a  uniform  time  step  At.  By  evaluating  the  cross-deriva¬ 
tive  terms  in  this  manner,  the  block  tridiagonal  form  of  the 
final  equations  is  maintained. 

Substituting  Eqs.  (50),  (52),  (54),  (57),  and  (59)  into 
Eq .  (48)  yields 


[/]  + 


0,  At 

1  +  0, 


+  |;(Ufl  -  IQ!  +[Sy}V-f>[S)n 


A”U 


At 


1  +  02 

,  fliAf 
1  +  02 

r, 

+  o 


^(-E  +  V,  +V2r +  ^(-F  +  W,  +Wa)« 


by 

(An~ 1 V2 )  +  ~  (An~  *  W, ) 


e2 


l  +e. 


An~  1  U 


(«i  (A/)a.(Af)9 


(61) 


where  [I]  is  the  identity  matrix.  In  Eq .  (61),  expressions 

such  as 

[-£(  [A]  -  [P]  *  [i?*]  )  n]  AnC/  (62) 

are  equivalent  to 

[-£(  U]  -  [i^]  +  )  "AnL73  (  63  ) 
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Step  3 : 


Unn  =  Un+  A  nU  (  68  ) 

In  Step  1,  An  U,  represents  the  remaining  terms  on  the  left- 
hand  side  of  Eq .  (64).  Equations  (66)  and  (67)  represent 

systems  of  equations  which  have  a  block  tridiagonal  structure. 
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For  the  two-dimensional  compressible  Navier-Stokes  equations 
the  block  matrices  have  a  dimension  of  4x4.  Central 
differencing  is  used  for  the  second  order  space  derivatives. 
Dissipation  terms  for  numerical  stability  are  added. 

Warming  and  Beam  [Ref.  15]  have  shown  that  the  algorithm 
can  be  simplified  by  assuming  that  p  is  locally  constant  so 
that  0p/0x=O,  0p/0y=O.  Then  [P]-[/?x]=0  and  [<?]-[  S  ]  =0 . 
Therefore,  the  algorithm  for  01  -hi,  0^  =  0,  implying  second 
order  accuracy  in  time,  obtains  the  following  form: 

[2]  +  •—  [-—  [A]  n~  [Rx3  n  A  n  Ux  =  RHS  (69) 

[I]  [-£;[B]  n--t[Syl  nAn  U  =  An  U,  (70) 

In  the  present  work  the  viscous  terms  were  treated  explicitly 
in  order  to  avoid  the  expensive  computation  of  the  matrices  Rx 
and  Sy,  The  Beam-Warming  algorithm  with  explicit  treatment  of 
viscous  terms  in  generalized  coordinates  ( I;  ,  )  is  written  as 

(  [I]  +—  [-?[An+J-Bn])AUn'1  = 

2  6i  *  (71) 

A  t{--!rFn--?rGn+4tRn+-zrS  ")  =  MS  n 

dr\  &r\ 

The  algorithm  is  further  simplified  by  approximately  factoriz¬ 
ing  the  LHS(61)  or  LHS(71)  operator  in  order  to  avoid  integra¬ 
tion  of  the  full  two-dimensional  operator.  The  factored  form 
of  the  algorithm  is 
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(  [I]  *-~dlAi°tk)  A  Ujft  =  RHS  (61)  n 


(72) 


( tj]  +^£6^"*)  au?;}  =  A  u;;i 

In  practice,  implicit  algorithms  have  stability  limits  due  to 
nonlinearities.  In  addition,  whenever  discrete  methods  are 
used  to  compute  high  Reynolds  number  viscous  behavior,  small 
scales  of  motion  appear  which  cannot  be  resolved  by  the 
numerics.  These  scales  are  brought  about  by  the  nonlinear 
interactions  in  the  convective  terms  of  the  momentum  equa¬ 
tions.  In  any  finite  discrete  mesh  the  small  scales  which 
cannot  be  resolved,  result  eventually  in  inaccuracy  and 
contamination  of  the  long  wavelength,  large  scale  phenomena. 
In  order  to  dampen  the  high  frequency  numerical  effects  caused 
by  the  poor  resolution  of  the  small  scales,  implicit  and 
explicit  numerical  dissipation  is  added  to  complete  the 
algorithm.  The  numerical  dissipation  terms  introduce  an  error 
level  that  does  not  interfere  with  the  accuracy  and  resolution 
of  any  physical  effects.  The  dissipation  terms  used  in  the 
present  work  are  the  ones  suggested  by  Ref.  15.  The  complete 
factorized  form  of  the  numerical  algorithms  with  the  implicit 
and  explicit  dissipation  terms  is 

U-  (it)  WtA”*(rJwI)e)]xtJ*  (it)  0„B»* 

-  At  (-3  ,F"-a,G”  +  3,R"  +  a,S”-e.wJD„) 
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C.  THE  BALDWIN-LOMAX  EDDY-VISCOSITY  MODEL 

The  Baldwin-Lomax  model  [Ref.  16]  is  an  algebraic  eddy 
viscosity  model  for  calculating  two-  and  three-dimensional 
separated  flows.  As  opposed  to  the  classical  boundary- layer 
approximation  which  assumes  zero  normal  pressure  gradient  in 
the  boundary  layer,  thereby  neglecting  the  normal  momentum 
equation,  the  Baldwin-Lomax  Thin  Layer  model  retains  the 
momentum  equations  and  makes  no  assumptions  about  the  pressure 
gradient.  The  advantage  of  this  model  arises  in  application 
to  high  Reynolds  number,  separated  turbulent  flows,  including 
reverse  flow  regions.  The  Baldwin-Lomax  model  is  a  two -layer 
algebraic  eddy-viscosity  model  and  avoids  the  difficulty  of 
finding  the  edge  of  the  boundary  layer.  The  effects  of 
turbulence  are  simulated  in  terms  of  an  eddy  viscosity 
coefficient  pt,  where  p+pt  replaces  p  (the  physical  viscosity 
coefficient)  in  the  stress  terms  of  the  governing  equations. 
In  the  model,  pt  is  given  by 

f  ^  inner  Y^Y c  .  . 

f*t  =  (74> 

l  t)  outer  YC>Y 

where  y  is  the  normal  distance  from  the  wall  and  yr  is  the 
smallest  y  where  inner  and  outer  values  are  equal.  In  the 
inner  region,  the  Prandtl-Van  Driest  formula  is  used: 

c>  inner  =  P  *  M  (75) 
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where  1  =  ky[l  -  exp(-y+/A  +  )]  and  |w|  is  the  magnitude  of  the 
vorticity.  In  two  dimensions, 


and 


M 


i  Ay  _  Ay  \ 

dy  dx 


(76) 


The  subscript  w  indicates  wall  values. 

For  the  outer  region 

outer  “  ^i?p  P  trr  vr>{y')  (78) 

where  K,  Ccp  are  constants  and 

F WAKE  =  ^MAX  fmax  for  boundary  layers,  or 

fWAK£  ”  cwk  y«AX  uZ  DIP/FMAX  for  wakes  and  separated  boundary 

layers. 

The  quantities  yMAX  and  FMAX  are  determined  from 

F(y)  =  y|v|  [1-exp  (y4/A  +  )  ] 

The  exponential  term  in  this  equation  is  set  equal  to  zero  for 
wakes . 

Fmax  is  the  maximum  values  of  F(y),  which  occurs  at  a 
value  of  y  =  yMAX.  The  function  FKLEB(y)  is  the  Klebanoff 
intermittency  factor  given  by 
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The  quantity  £7qIF  is  the  difference  between  U  =  U  |  y=yMA)(  and 
minimum  total  velocity  at  a  fixed  flow-wise  station: 


(u2  +  v2 


ru2+  v2  ). 


(81) 


where  the  second  term  is  zero  except  in  wakes. 

In  order  to  achieve  agreement  with  Cebeci  [Ref.  17],  the 
values  determined  for  the  constants  are: 


A + 

=  26 

k 

=  0.4 

C,P 

=  1.6 

K 

=  0.0168 

C  KLEB 

=  0.3 

Pr 

=  0.72 

C  WK 

=  1.0 

pr  t 

=  0.9 

D  BOUNDARY  CONDITIONS 

The  following  boundary  conditions  were  used  in  the 
numerical  implementation.  A  non-slip  non-penetration  boundary 
condition  in  terms  of  the  contravar iant  velocity  components 
was  used  on  the  airfoil  solid  boundary.  The  unsteady  pitching 
motion  was  accomplished  by  rotating  the  grid  about  the  quarter 
chord  point  at  pitch  rates  determined  by  the  desired  reduced 
frequency.  The  inflow  boundary  was  placed  approximately  eight 
chord  lengths  away  from  the  body  surface.  Therefore, 
freestream  conditions  were  specified  at  the  grid  inflow 
boundary.  Simple  first  order  extrapolation  was  used  for  the 


flow  variables  at 

the 

outflow  boundary. 

For 

the 

wake , 

averaging  of  the 

f  low 

variables  on  the 

upper 

and 

lower 

surfaces  of  the  C-grid  was  used.  The  velocity  on  the  surface 
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is  given  by  us=x=ti)2  and  vs=z=-ox.  The  contravar iant  velocity 
components  for  viscous  flow  solutions  are  set  equal  to  0. 
Therefore,  the  non-slip  condition  in  terms  of  the  physical 
velocity  components  for  moving  grid  is  expressed  by 


u 

1 

r\y  -5 , 

X 

V. 

J 

±. 

E.  GRID  GENERATION 

1.  Grid  Generation  Methods 

In  order  to  effectively  compute  complex  flowfields,  the 
physical  domain  of  interest  must  be  discretized  with  a  finite 
mesh.  The  requirements  of  an  efficient  numerical  grid  are 
[Ref.  18]: 

1.  smooth  grid  lines  so  that  the  transformation  derivatives 
(metrics)  are  continuous. 

2.  grid  point  spacing  which  varies  inversely  with  expecta¬ 
tion  of  large  numerical  errors. 

3.  minimizing  grid  skewness  to  avoid  large  truncation 
errors . 

Several  general  grid  generation  techniques  exist;  among  these 
the  most  common  methods  are: 

1.  Complex  Variable  methods,  where  the  transformations  are 
at  least  partly  analytic;  this  method  is  restricted  to 
two- dimensions . 

2.  Algebraic  methods,  usable  in  two-  or  three-dimensions. 

3.  Differential  Equation  methods,  usable  in  two-  or  three- 
dimensions  . 
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The  grid  generation  method  used  in  this  study  utilized  the 
algebraic  grid  generation  technique  because  of  the  computa¬ 
tional  efficiency  and  spe_d  it  provides. 

2.  Algebraic  Grid  Generation  Method 

The  algebraic  method  used  employs  known,  easily 
invertible,  functions  to  map  arbitrarily  shaped  regions  (in 
this  case,  the  airfoil  contour)  into  a  simpler  computational 
domain.  The  airfoil  surface  is  unwrapped  to  form  a  simple 
curve  in  the  computational  plane  as  in  Figure  6.  In  the 
computational  domain  lines  are  first  drawn  normal  to  this 
curve,  and  the  grid  points  in  the  normal  direction  are 
subsequently  generated.  In  the  case  of  airfoil  grid  genera¬ 
tion,  consideration  must  be  given  for  the  clustering  of  grid 
points  near  the  airfoil  in  order  to  adequately  resolve  the 
near  surface  viscous  layers.  An  algebraic  function  is  used  to 
provide  a  uniform  stretching  normal  to  the  airfoil  in  the 
computational  plane.  The  resulting  computational  domain  grid 
is  wrapped  back  to  the  initial  physical  grid  using  inverse 
transformations . 

The  grid  generation  program  used  in  this  study  yielded 
a  157x58  C-type  grid.  The  program  is  listed  in  Appendix  C. 
Grids  for  the  airfoils  under  consideration  are  displayed  in 
Figures  7-9.  Figure  7  displays  the  global  grid  including  wake 
for  the  NACA  0012  airfoil.  Figures  8  and  9  show  the  body- 
fitted  grid  in  more  detail  for  the  two  modified  airfoils. 
Grid  clustering  near  the  airfoils  is  shown  for  resolving  the 
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Figure  6.  Airfoil  Grid  Unwrapping 


boundary  layer  flow  in  turbulent  and  separated  flow  condi¬ 
tions  . 
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Figure  7.  N0012  Global  Grid 
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Figure  8.  N0012-63  Local  Grid 


Figure  9.  N0012-33  Local  Grid 
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IV.  RESULTS  AND  DISCUSSION 


A.  INVESTIGATION  METHOD 

Two-dimensional  unsteady  flows  involving  both  harmonically 
oscillating  and  rapidly  pitching  (ramp  pitch-up)  airfoils  were 
investigated.  Studies  for  the  oscillatory  cases  were  conduct¬ 
ed  using  the  NACA  0012  airfoil  in  order  to  enable  comparisons 
with  previously  conducted  experimental  measurements.  Flows 
over  all  three  airfoils  rapidly  pitched  from  0°  to  30°  angle 
of  attack  with  the  reduced  frequencies  and  Mach  numbers  (shown 
in  Table  2)  were  subsequently  investigated. 

TABLE  2.  RAPIDLY  PITCHING  AIRFOIL  CASES,  Re  =  4  X  10 6 


AIRFOIL 

0.3 
k= .  01 

Mach 
k= .  02 

0 . 4  Mach 
k= . 01  k= . 02 

NACA 

0012 

X 

X 

X 

X 

NACA 

0012-63 

X 

X 

X 

X 

NACA 

0012-33 

X 

X 

X 

X 

Solutions  for  the  harmonically  oscillating  NACA  0012 
airfoil  were  obtained  for  flows  at  a  freestream  Mach  number  of 
0.3,  Reynolds  number  based  on  the  root  chord  Rec=4xl06  and  for 
two  reduced  frequencies  of  k=0 . 1  and  k=0 . 2 .  The  reduced 
frequency  is  defined  as  k=o)c/2Us  for  the  oscillatory  case 
where  a(t)=a0  +  a.  sin(wt) .  In  terms  of  nondimensional  quanti¬ 
ties  k  is  given  by  k=(i)/2Ma,  and  &(t)=a,  cocosfot)  is  the 
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instantaneous  pitch-up  rate  which  varies  during  the  cycle  with 
0)-2kM#|  for  a  unit  chord  length.  The  present  computations  were 
conducted  for  a  variation  of  the  angle  of  attack  as  «(t)=10°+ 
6sin(t).  Experimental  test  conditions  for  the  same  freestream 
Mach  and  Reynolds  numbers  were  for  a(t)=10°+  5sin(t) .  Because 
the  airfoil  stalled  just  at  15°  for  the  experiment,  the 
computed  flow  for  the  same  conditions  did  not  yield  a 
hysteresis  loop.  Therefore,  as  McCroskey  suggested,  the 
oscillation  was  increased  by  one  degree  and  a  hysteresis  loop 
was  obtained.  The  computed  lift  behavior  shown  in  Figures  10 
and  11  exhibits  the  well-known  hysteresis  loop  of  a  harmoni¬ 
cally  oscillating  airfoil  experiencing  dynamic  stall.  The 
computation  was  initiated  from  a  steady-state  solution  at  a=5° 
and  was  carried  out  for  two  cycles.  Figures  10  and  11  display 
the  results  of  the  second  cycle  from  the  two-cycle  computa¬ 
tion. 

Flow  solutions  for  a  rapidly  pitching  airfoil  were 
obtained  by  pitching  the  airfoil  at  a  constant  rate  from  a 
zero  angle  of  attack  and  steady-state  flow  conditions  to  an 
angle  of  attack  of  30°  at  the  desired  reduced  frequency  and 
freestream  Mach  nu;  ?r.  For  the  case  where  ramp  motion  was 
imposed,  the  reduced  frequency  k  is  given  by  k=&c/2Uw  where  a 
is  the  constant  pitch-up  rate.  In  terms  of  nondimensional 
quantities  k=<i>/2Ma,  or  w=2kMa)=constant ,  and  a  ( t )  =a  5+  (  a  -  a  c)  ut , 
where  in  this  study  ag=08  and  ai=30°  and  <i>=&(t).  A  summary  of 
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the  computed  results  for  the  rapidly  pitching  airfoils  is 
provided  in  Table  3. 

TABLE  3.  PEAK  LIFT  COEFFICIENTS 


AIRFOIL 

Mach 

No  . 

Reduced 

Frequency 

Peak 

C1 

Angle  of 
Attack 

0 . 3 

.01 

1 . 90 

19.40° 

N0012 

.  02 

2  .  10 

23.40° 

0 . 4 

.01 

1 .97 

21.08° 

.02 

2 . 24 

25.21° 

0 . 3 

.01 

1 . 84 

18.50° 

N001 2-63 

.02 

2 . 09 

23.00° 

0 . 4 

.01 

1 . 93 

20.60° 

.  02 

2 . 22 

25.00° 

0 . 3 

.01 

1 . 65 

15.80° 

N0012-33 

.02 

1 . 94 

19.60" 

0.4 

.01 

1 . 74 

17.10° 

.  02 

2 . 09 

23.15° 

A  general  observation  on  the  computed  solution  is  that  the 
modified  NACA  0012-63  airfoil  exhibited  comparable  lift 
behavior  to  the  baseline  NACA  0012  airfoil.  A  slight  angle  of 
attack  versus  lift  curve  shift  was  observed  between  the  two 
airfoils  at  all  reduced  f requency/Mach  number  combinations. 
However,  for  the  same  flow  parameters  the  NACA  0012-33  airfoil 
consistently  underperformed  the  larger  leading  edge  radius 
airfoils.  The  difference  in  the  unsteady  lift  behavior  with 
increasing  angle  of  attack  resulted  from  the  different  flow 
character  at  the  leading  edge  region  as  will  be  shown  later  in 
detail.  In  the  following  sections,  detailed  comparison  of  the 
lift  behavior  of  the  three  airfoils  'is  presented  at  the 
various  flow  conditions.  The  flow  characteristics  and  the 
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development  and  progression  of  the  dynamic  stall  process  are 
examined.  The  effects  of  reduced  frequency  and  freestream 
Mach  number  are  also  discussed. 

B.  LIFT  BEHAVIOR 

1.  Harmonically  Oscillating  Airfoil 

Although  hysteresis  loops  characteristic  ».  c  an 
oscillating  airfoil  undergoing  dynamic  stall  were  observed  for 
reduced  frequencies  of  both  0.1  and  0.2,  agreement  with 
experimental  results  varied.  At  a  reduced  frequency  of  0.1, 
the  response  agreed  well  with  McCroskey's  experimental  results 
(Figure  10)  during  the  upstroke.  Maximum  lift  coefficient  was 
1.52  at  an  angle  of  attack  of  15.8°.  The  lift  behavior 
continues  to  agree  well  with  the  experimental  data  during  the 
upstroke  of  the  hysteresis  loop  and  during  the  initial  part  of 
the  downstroke.  As  the  downstroke  continues  and  the  flow 
reattaches,  however,  the  numerical  solution  displays  signifi¬ 
cantly  greater  oscillations  in  lift  coefficient  compared  to 
experimental  data,  which  were  obtained  as  an  average  over 
several  cycles  of  oscillation. 

At  a  reduced  frequency  of  0.2  (Figure  11),  the  numeri¬ 
cal  solution  exhibited  a  much  smaller  hysteresis  loop  than  the 
respective  experimental  data.  In  this  case,  maximum  lift 
coefficient  in  the  numerical  solution  occurred  just  prior  to 
the  downstroke,  so  that  flow  reattachment  occurred  almost 
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Figure  10.  Lift  Behavior,  NACA  0012,  M=0 . 3 


immediately.  Maximum  lift  coefficient  at  k=0.2  was  only 
slightly  higher  than  at  the  k=0.1  case. 

The  lack  of  agreement  with  experimental  data  at  a 
reduced  frequency  of  0.2  and  during  the  flow  reattachment 
process  at  k=0.1  may  be  indicative  of  the  poor  behavior  of  the 
eddy-viscosity  model  (which  is  suitable  for  steady  flows). 
Higher  values  of  the  reduced  frequency  resulted  in  larger 
discrepancies  from  the  measured  lift  values. 

2.  Rapidly  Pitching  Airfoil 

Lift  coefficient  vs.  angle  of  attack  comparison  of  the 
three  airfoils  is  presented  in  Figures  12-15  for  the  Mach 
number /reduced  frequency  combinations  listed  in  Table  2.  For 
the  same  flow  parameters,  all  three  airfoils  have  nearly  the 
same  lift  curve  slope  until  the  onset  of  dynamic  stall.  Along 
with  the  aforementioned  angle  of  attack  shift  of  the  lift 
curves  between  the  NACA  0012  and  -63  airfoils  after  stall 
occurs,  the  peak  lift  coefficients  for  the  NACA  0012  airfoil 
are  slightly  higher  than  for  the  -63  airfoil.  The  peak  lift 
coefficient  for  the  NACA  0012  occurred  0.2  to  0.9  degrees 
angle  of  attack  higher  than  the  peak  CL  for  the  -63  airfoil 
for  the  same  flow  conditions.  With  both  airfoils  having  the 
same  leading  edge  radius,  the  small  difference  in  performance 
may  be  attributed  to  the  slightly  thicker  contouring  of  the 
forward  part  of  the  baseline  NACA  0012  airfoil  (Table  1). 
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Figure  12.  Lift  Comparison,  M=0.3,  k=0.01 
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Ramp  pitch-up  from  0  to  30  deg,  M=0.3,  Re=4*  10**6, 


Figure  13.  Lift  Comparison,  M=0.3,  k=CK02 


48 


Ramp  pitch-up  from  0  to  30  deg,  M=0.4,  Re=4*  10**6, 


Figure  14.  Lift  Comparison,  M=0.4,  k=0.01 
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Ramp  pitch-up  from  0  to  30  deg,  M=0.4,  Re=4*  10**6, 


Figure  15.  Lift  Comparison,  M=0 . 4 ,  k=0.02 
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The  maximum  lift  coefficient  obtained  with  the  NACA 


0012-33  airfoil  was  lower  than  that  obtained  with  the  other 
two  airfoils  and  occurred  at  1.8  to  4.0  degrees  lower  angle  of 
attack.  At  a  reduced  frequency  of  k=0.01  and  freestream  Mach 
number  of  0.3,  the  initial  leading  edge  vortex  originated  at 
an  angle  of  attack  of  12°  for  the  -33  airfoil.  Under  the  same 
flow  conditions,  formation  of  the  vortex  occurred  at  17°  and 
.16.5°  for  the  baseline  0012  and  -63  airfoils,  respectively. 
Similar  differences  occurred  at  all  other  flow  conditions 
investigated.  The  smaller  leading  edge  radius  of  the  NACA 
0012-33  airfoil  promoted  earlier  development  of  the  dynamic 
stall  vortex  from  the  leading  edge  upper  surface. 

As  will  be  shown,  the  dynamic  stall  vortex  originates 
as  a  result  of  the  combination  of  the  accelerated  flow  over 
the  leading  edge  with  the  boundary  layer  reverse  flow  behind 
the  leading  edge.  An  adverse  pressure  gradient  is  encountered 
by  the  flow  as  it  passes  the  suction  peak  just  downstream  of 
the  airfoil  leading  edge.  This  adverse  pressure  gradient 
causes  flow  deceleration  just  aft  of  the  suction  peak,  leading 
eventually  to  the  boundary  layer  separating  at  a  critical 
value  of  the  adverse  pressure  gradient.  The  momentum  of  the 
flow,  however,  is  increased  as  the  flow  accelerates  around  the 
leading  edge  and  opposes  the  tendency  of  the  flow  to  separate. 
Eventually,  depending  on  the  airfoil  shape  characteristics, 
the  boundary  layer  separates,  reversed  flow  occurs,  and 
combines  with  the  freestream  to  form  thfc  leading  edge  vortex. 
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As  observed  in  this  investigation,  the  size  of  the 
leading  edge  radius  is  of  primary  importance  in  determining 
the  angle  of  attack  at  which  the  boundary  layer  separates  and 
rolls  to  form  the  dynamic  stall  vortex.  For  the  same  flow 
parameters  and  angle  of  attack  during  pitch-up,  the  critical 
pressure  gradient  for  flow  separation  occurred  earlier  on  the 
airfoil  with  smaller  leading  edge  radius  (NACA  0012-33). 
Stated  otherwise,  at  a  given  angle  of  attack  during  pitch-up, 
the  adverse  pressure  gradient  aft  of  the  suction  peak  is 
greater  for  the  smaller  leading  edge  airfoil. 

The  contouring  of  the  forward  part  of  the  airfoil  is 
of  secondary  importance  in  alleviating  development  of  the 
adverse  pressure  gradient.  The  effect  of  contouring  is 
indicated  by  the  slightly  earlier  development  of  leading  edge 
reverse  flow  on  the  NACA  0012-63  airfoil  compared  to  the 
baseline  0012  airfoil. 

In  summary,  enlarging  the  leading  edge  radius  and 
thickening  the  contouring  of  the  forward  part  of  the  airfoil 
results  in  decreasing  the  adverse  pressure  gradient  encoun¬ 
tered  by  the  flow,  with  resulting  delay  in  separation.  This 
is  shown  graphically  in  Figures  16-20  where  instantaneous 
streamlines  are  shown  for  the  three  airfoils  at  the  same  17° 
angle  of  attack,  0.4  Mach,  and  k=0.01.  Reversed  boundary 
layer  flow  has  just  begun  on  the  NACA  0012  airfoil  aft  of  the 
suction  peak,  whereas  initiation  of  the  dynamic  stall  vortex 
has  already  occurred  on  the  NACA  0012-63  airfoil  which  has 
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thinner  contouring.  The  NACA  0012-33  airfoil  with  the  small 
leading  edge  radius  has  a  fully  developed  dynamic  stall  vortex 
and  is  0.6°  angle  of  attack  prior  to  dynamic  stall.  A 
detailed  discussion  of  the  vortical  flow  development  follows. 
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C.  VORTEX  FLOW  DEVELOPMENT 


1.  Rapidly  Pitching  Airfoil 

The  development  of  the  vortical  flowfield  as  the  angle 
of  attack  is  increasing  follows  a  consistent,  general  sequence 
of  events.  This  sequence  of  events  requires  examination  of 
the  variation  of  several  flow  quantities  as  angle  of  attack  is 
increased.  Figures  21-33  illustrate  the  flow  development 
stages  using  the  velocity  vector  field,  the  pressure,  and 
vorticity  field  contours  for  the  NACA  0012-63  airfoil  at  M=0.4 
and  k=0.01.  Figures  34-47  display  the  associated  surface 
pressures  and  temperatures  of  the  developing  flow. 

Initially,  smooth  streamlined  flow  is  observed  over  the 
airfoil  leading  edge  (Figures  21-22).  As  angle  of  attack 
increases  beyond  a  certain  critical  limit  depending  on 
pitching  rate,  freestream  Mach  number,  airfoil  shape,  and 
Reynolds  number,  small  reverse  flow  regions  develop  on  the 
upper  surface  aft  of  the  suction  peak  and  forward  of  the 
trai} ing  edge  due  to  the  adverse  pressure  gradients  encoun¬ 
tered  by  the  flow.  Figure  23  shows  the  reversed  boundary 
layer  flow  near  the  airfoil  surface.  Figure  24  displays  the 
developing  leading  edge  vortex. 

The  developing  vortex  forms  on  the  upper  surface  just 
aft  of  the  leading  edge  as  a  result  of  the  combination  of  the 
accelerated  flow  near  the  suction  peak  and  the  reverse  flow 
j  >  t  aft  of  the  suction  peak.  It  is  observed  that  the  dynamic 
stall  vortex  initiates  as  a  separated  flow  bubble,  and  rapidly 
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grows  in  size  to  form  the  characteristic  leading  edge  dynamic 
stall  vortical  structure  as  the  angle  of  attack  increases. 
During  this  initiation  process  the  primary  vortex  has  a  oval 
shape  and  becomes  more  rounded  as  the  angle  of  attack  increas¬ 
es.  Figure  25  shows  the  oval  dynamic  stall  vortex  taking 
shape.  Figure  26  displays  the  relative  strength  of  the 
developing  vortex. 

The  vortex  grows  in  size  and  strengthens,  becoming 
approximately  circular  as  angle  of  attack  increases,  and  moves 
away  from  the  airfoil.  Figure  27  shows  the  development  of  the 
vortex  as  it  moves  downstream.  Passage  of  the  dynamic  stall 
vortex  over  the  airfoil  surface  induces  reverse  velocities  and 
significantly  contributes  to  the  development  of  reversed 
flows.  A  secondary  vortex  originates  near  the  leading  edge 
and  also  grows  as  angle  of  attack  increases.  Figures  28  and 
29  show  the  development  of  a  secondary  vortex  as  the  primary 
vortex  moves  downstream  and  promotes  reverse  flow  over  the 
whole  upper  airfoil  surface. 

At  the  trailing  edge,  flow  separation  is  initiated 
at  approximately  the  same  time  as  at  the  leading  edge,  and  a 
pressure  gradient  exists  between  the  lower  and  upper  surfaces. 
The  upper  surface  reverse  flow  combines  with  the  high  pressure 
flow  from  the  lower  surface  to  form  a  counter-clockwise  vortex 
at  the  trailing  edge  shortly  before  dynamic  stall  occurs. 
Figures  30  and  31  show  the  developed  trailing  edge  vortex  0.6° 
angle  of  attack  before  peak  lift  coefficient  is  obtained.  As 
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the  primary  vortex  grows  and  moves  downstream  past  approxi¬ 
mately  60-70%  chord,  dynamic  stall  occurs,  usually  accompanied 
by  development  of  a  small,  tertiary  vortex  at  the  leading 
edge.  Figures  32  and  33  display  the  position  of  the  dynamic 
stall  vortex  0.4°  angle  of  attack  after  attainment  of  peak 
lift  coefficient. 

By  this  time,  the  primary  vortex  is  centered  at  a 
distance  greater  than  the  airfoil  maximum  thickness  from  the 
airfoil  and  combines  with  lower  surface  high  pressure  flow  to 
energize  the  counter-clockwise  vortex  at  the  trailing  edge. 

With  further  increase  in  angle  of  attack,  the  airfoil 
enters  the  deep  stall  regime,  and  the  freestream  flow  has  no 
effect  on  the  upper  surface  flow  characteristics.  In  this 
case,  the  aerodynamic  behavior  of  the  airfoil  is  greatly 
determined  by  the  vortical  flow  field  formed  by  the  shed 
vortices,  and  the  strength  of  the  vortices  themselves  deter¬ 
mine  the  flow  over  the  upper  surface.  Figures  48  and  49  show 
the  NACA  0012-33  airfoil  in  deep  stall  at  0.3  Mach,  k=0.02, 
a=22.0B  and  demonstrate  that  the  vortices  can  combine  to  form 
a  large  vortical  region  and  entrain  other  secondary  flows.  In 
this  case  a  counterclockwise  vortex  was  formed  between  the 
paired  primary  and  secondary  vortices  and  the  airfoil. 
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Figure  23.  Velocity  Field,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
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Figure  25.  Velocity  Field,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
a=17 . 0® 
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Figure  27.  Velocity  Field,  NACA  0012-63,  M=0.4,  k=0.01, 
a=18 . 0° 
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Figure  28. 


Velocity  Field,  NACA  0012-63,  M=0 . 4 .  k=0.01, 
a  =  19 . 0® 
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Figure  34.  Surface  Pressure,  NACA  0012-63,  M=0.4,  k=0.01, 
a=14 . 0° 
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Chordwise  Position,  xJc 
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Figure  35.  Surface  Pressure,  NACA  0012-63,  M=0.4,  k= 
a  =  l 6 . 0 0 
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Chordwise  Position,  x/c 


Figure  37.  Surface  Pressure,  NACA  0012-63,  M=0.4,  k=0.01, 
et  =  l  8 . 0  0 
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Chordwise  Position,  x/c 


N0012-63  Airfoil,  M=0.4,  Re=4*10* 


Figure  38.  Surface  Pressure,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
«=19.0° 
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N0012-63  Airfoil,  M=0.4,  Re=4*10**6,  k=0.01,  a=20.0  deg 


Figure  39.  Surface  Pressure,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
a  =  20 . 0 0 
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N0012-63  Airfoil,  M=0.4,  Re=4*10**6,  k=0.01,  a=21.0  deg 


Figure  40.  Surface  Pressure,  NACA  0012-63,  M=0.4,  k=0.01, 
«=21.0° 
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NOO 12-63  Airfoil,  M=0.4,  Re=4*10**6,  k=0.01,  a=14.0  deg 


Figure  41.  Surface  Temperature,  NACA  0012-63,  M=0.4,  k=0.01, 
a=14 . 0° 
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2-63  Airfoil,  M=0.4,  Re=4*: 
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Figure  42.  Surface  Temperature,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
a  =  1 6 . 0  9 
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Figure  43.  Surface  Temperature,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
a  =  l 7 . 0 0 
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NOO 12-63  Airfoil,  M=0.4,  Re=4*10**6,  k=0.01,  a=18.0  deg 


Figure  44.  Surface  Temperature,  NACA  0012-63,  M=0.4,  k=0.01, 
<x  =  18 . 0° 
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=0.4,  Re=4*  10**6,  k=0.01,  a=19.0  deg 
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Figure  45.  Surface  Temperature,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
a=19.0° 
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=0.4,  Re=4*  10**6,  k=0.01,  a=20.0  deg 
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Figure  46.  Surface  Temperature,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
o=20.0° 
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Chordwise  Position,  x/c 


N0012-63  Airfoil,  M=0.4,  Re=4*10**6,  k=0.01,  a=21.0  deg 
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Figure  47.  Surface  Temperature,  NACA  0012-63,  M=0 . 4 ,  k=0.01, 
0=21.0° 
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Chordwise  Position,  x/c 
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Figure  < 


2.  Harmonically  Oscillating  Airfoil 

The  process  of  vortical  flow  development  at  the  leading 
edge  during  the  upstroke  is  similar  to  that  of  a  rapidly 
pitching  airfoil  examined  in  the  previous  section.  Although 
the  pitching  rate  is  not  constant  in  the  case  of  the  harmoni¬ 
cally  oscillating  airfoil,  boundary  layer  separation  and 
vortex  development  and  shedding  also  occur  during  the  up¬ 
stroke  . 

The  flow  reattachment  process  during  the  downstroke 
involves  the  shedding  of  the  large  primary  vortex  into  the 
wake  and  the  diminishing  intensity  of  the  trailing  edge 
vortex.  Prior  to  the  downstroke,  the  combination  of  the  large 
clockwise  primary  vortex  and  the  counter-clockvdse  trailing 
edge  vortex  cause  extensive  reverse  flow,  even  outside  the 
boundary  layer,  over  the  airfoil  upper  surface  aft  of  30% 
chord  (Figure  50).  As  the  angle  of  attack  decreases  further, 
flow  over  the  leading  edge  upper  surface  rapidly  reattaches, 
while  the  primary  vortex  is  centered  above  the  trailing  edge 
(Figure  51).  With  further  reduction  in  angle  of  attack,  the 
primary  vortex  moves  downstream  of  the  airfoil,  and  the 
trailing  edge  vortex  diminishes  (Figure  52).  Reverse  flow  at 
this  time  exists  only  on  aft  portions  of  the  airfoil.  As  the 
angle  of  attack  approaches  10-11°,  the  primary  vortex  has  been 
swept  downstream,  the  trailing  edge  vortex  has  diminished,  and 
smooth,  attached  flow  is  established  over  the  upper  surface 
[ Figure  53 ) . 
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Figure  50.  Velocity  Field,  NACA  0012,  M=0.3,  k=0.1,  r= 
downstroke 
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Figure  52.  Velocity  Field,  NACA  0012,  M=0 . 3 ,  k=0 . 1 ,  a=11.0 
downstroke 
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Figure  53.  Velocity  Field,  NACA  0012,  M=0 . 3 ,  k=0 . 1 
downstroke 


D.  REDUCED  FREQUENCY  EFFECT 

Al  a  given  Mach  number  the  three  airfoils  exhibited  higher 
peak  lifr  coefficients  at  respectively  higher  angles  of  attack 
at  a  reduced  frequency  of  k=0.02  when  compared  to  tnat  of 
k=0.01.  Lift  curve  slopes  for  a  given  airfoil  at  the  two 
reduced  frequencies  were  nearly  identical  to  the  point  of 
initia'  flow  separation  and  subsequent  dynamic  stall  at 
k=0.01.  Lift  coefficient  at  k=0.02  continued  increasing, 
showing  the  increase  in  unsteady  lift  with  higher  pitching 
rate.  The  effect  of  increasing  the  reduced  frequency  at  a 
freestream  Mach  number  of  0.3  is  shown  in  Figures  54  and  55 
for  the  NACA  0012  and  NACA  0012-33  airfoils,  respectively.  At 
the  higher  reduced  frequency  of  0.02,  formation  of  the  primary- 
vortex  occurred  at  angles  of  attack  0.7°  to  1.6°  higher  than 
at  a  reduced  frequency  of  0.01  in  all  cases.  As  a  result,  the 
entire  vortical  flow  development  was  delayed  to  progressively 
higher  angles  of  attack.  Dynamic  stalling  angle  occurred  4.0° 
to  4.8°  higher  at  the  higher  reduced  frequency  and  the 
resulting  peak  lift  coefficients  were  15-20%  higher  compared 
to  the  k=0.01  cases  (Table  3). 

Experimental  work  by  Chandrasekhara  and  Carr  [Ref.  3] 
using  the  NACA  0012  airfoil  demonstrated  the  same  trends, 
wherein  a  higher  reduced  frequency  resulted  in  retaining  the 
vortex  above  the  airfoil  to  higher  angles  of  attack,  thereby 
delaying  dynamic  stall. 
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Figure  54.  Reduced  Frequency  Effect,  NACA  0012,  M=0.3 
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Ramp  pitch-up,  N0012-32 


Figure  55.  Reduced  Frequency  Effect,  NACA  0012-33,  M=0 . 3 
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The  effect  of  the  increased  pitching  rate  in  delaying  the 
onset  of  dynamic  stall  and  increasing  the  unsteady  1 ifi  may 
prove  to  be  of  practical  operational  benefit  as  knowledge  of 
the  details  of  the  flow  development  increases. 

E.  EFFECT  OF  FREESTREAM  MACH  NUMBER 

The  effect  of  increasing  freestream  Mach  number  from  0.3 
to  0.4  resulted  in  slightly  displacing  the  lift  curve  upward 
for  the  higher  Mach  number.  The  slope  of  the  lift  curve  was 
unchanged.  However,  for  a  given  angle  of  attack,  a  very 
s]jght  increase  in  lift  coefficient  was  observed  at  the  higher 
Mach  number  before  the  onset  of  dynamic  stall.  Figures  56  and 
57  show  the  effect  of  increasing  Mach  number  for  the  NACA  0012 
and  0012-33  airfoils  at  a  reduced  frequency  of  k-0.02.  Flow 
at  a  Mach  number  of  0.4  resulted  in  slightly  higher  peak  lift 
coefficients  than  at  M=0 . 3  .  In  addition,  the  peak  lift 
coefficient  at  M=0 . 4  occurred  1-2°  higher  angle  of  attack  than 
at  M=0.3. 

The  initial  development  of  reverse  flow  regions  and 
vorticity  was  largely  independent  of  Mach  number'  for  a  given 
airfoil  in  that  the  primary  vortex  originated  at  approximately 
the  same  angle  of  attack  for  the  two  Mach  numbers.  Subsequent 
growth  of  the  primary  vortex  and  development  of  the  secondary, 
tertiary,  and  trailing  edge  vortices  are  delayed  to  slightly 
higher  angles  of  attack  at  M=0 . 4  compared  to  the  M=0 . 3  case. 
As  a  result,  dynamic  stall  occurs  at  a  higher  angle  of  attack. 
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Ramp  pitch-up,  N0012  airfoil,  Re  4*  10**6,  k  =  0.02 


Ramp  pitch-up,  N0012-33  airfoil,  Re  4*  10**6,  k  =  0.02 


Ekaterinaris  [Ref.  11]  found  that  for  the  SSC-a90  airfoil 
with  9%  thickness  ratio,  increasing  Mach  number  from  0.2  to 
0.4  resulted  in  decreasing  both  lift  coefficient  and  dynamic 
stalling  angle  of  attack  at  a  reduced  frequency  of  0.01.  The 
experimental  work  by  Chandrasekhara  and  Carr  [Ref.  3]  using  a 
harmonically  oscillating  NACA  0012  airfoil  also  showed  that 
increasing  Mach  number  reduces  the  angle  of  attack  at  which 
dynamic  stall  occurred.  Further  experimental  and  computation¬ 
al  investigation  is  required  at  comparable  pitching  rates  to 
ascertain  the  accuracy  of  the  computational  model  used  at 
compressible  flow  Mach  numbers.  Further  investigatioj}  into 
the  influence  of  airfoil  thickness  ratio  on  Mach  number 
effects  is  needed. 

F.  PRESSURE  GRADIENT  AT  FLOW  SEPARATION 

The  peak  pressure  gradient  observed  at  initial  flow 
separation  occurring  aft  of  the  leading  edge  was  investigated 
for  the  three  airfoils  at  the  reduced  frequencies  and  Mach 
numbers  listed  in  Table  2.  Figure  58  displays  a  plot  of 
streamwise  pressure  gradient  along  the  airfoil  chord  for  the 
NACA  0012-63  airfoil  at  different  freestream  conditions.  The 
figure  shows  that,  for  the  same  airfoil,  the  peak  pressure 
gradient  encountered  by  the  flow  at  initial  flow  separation  is 
independent  of  freestream  speed  or  pitching  rate.  Similar 
results  were  obtained  for  the  other  airfoils.  Figure  59  shows 
the  streamwise  pressure  gradient  for  the  NACA  0012-33  airfoil. 
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Figure  60  displays  pressure  gradients  at  initial  flow 
separation  angle  of  attack  for  the  three  airfoils  at  the 
freestream  condition  of  0.4  Mach  and  reduced  frequency  of 
k=0.01.  The  peak  pressure  gradient  observed  at  the  instant  of 
flow  reversal  was  higher  for  the  NACA  0012-33  airfoil  than  for 
the  airfoils  with  larger  leading  edge  radius.  Similar  results 
were  obtained  at  other  freestream  conditions.  The  streamwise 
location  of  the  peak  pressure  gradient  for  the  NACA  0012-33 
airfoil  was  slightly  upstream  of  that  for  the  NACA  0012-63 
airfoil  as  shown  in  Figure  61.  The  higher  peak  pressure 
gradient  observed  for  the  NACA  0012-33  airfoil  is  an  indica¬ 
tion  that  flow  momentum  near  the  surface  is  higher  around  the 
smaller  leading  edge  radius.  The  higher  flow  momentum 
combined  with  the  relatively  upstream  location  of  the  peak 
pressure  gradient  and  sharper  flow  turning  angle  combined  to 
cause  the  magnitude  of  the  peak  pressure  gradient  for  the  NACA 
001  2-33  airfoil  to  be  33%  greater  than  that  of  the  NACA 
0012-63  airfoil.  Thus,  the  critical  pressure  gradient 
required  for  flow  reversal  is  dependent  on  leading  edge 
radius,  in  that  a  smaller  leading  edge  radius  increases  the 
peak  pressure  gradient  at  flow  reversal. 

The  streamwise  location  of  flow  separation  is,  in  ti’rn, 
dependent  on  the  location  of  the  peak  pressure  gradient.  The 
actual  location  of  flow  reversal  on  the  airfoil  surface 
occurred  1.0-1. 5%  chord  length  dov/nstream  of  the  location  of 
the  peak  pressure  gradient  as  shown  in  Figure  61 .  This  delay 
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is  attributed  to  the  time  lag  in  the  aerodynamic  response  to 
the  rapid  increase  in  pressure  gradient.  The  more  downstream 
location  of  the  peak  pressure  gradient  of  the  NACA  0012-63 
airfoil  caused  a  similar  relative  location  of  the  position  of 
initial  flow  reversal  compared  the  NACA  0012-33  airfoil. 


M=0.4,  k=0.01,  Re=4*10**6 
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Figure  60.  Pressure  Gradient  Comparison,  M=0 . 4 ,  k=0.01 
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Figure  61.  Flow  Reversal  Location,  M=0 . 4 ,  k=0.01 
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CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

1.  The  results  of  the  unsteady  flow  solutions  displayed  good 
agreement  with  the  exper imental  re  ults  for  the  harmoni¬ 
cally  oscillating  airfoil  at  a  reduced  frequency  of  0.1 
and  Mach  number  of  0.3.  At  a  higher  reduced  frequency  of 
0.2,  however,  "he  agreement  with  experimental  results  was 
poor  during  the  downstroke. 

2.  The  leading  edge  geometry  of  an  airfoil  was  found  to  have 
a  significant  effect  on  the  development  of  the  vortical 
flowfield  and  dynamic  stall  characteristics  of  the 
airfoil.  Of  primary  importance  is  the  size  of  the 
leading  edge  radius.  A  larger  leading  edge  radius  delays 
development  of  the  adverse  pressure  gradient  necessary 
for  boundary  layer  separation  and  eventual  vortex  forma¬ 
tion.  Of  secondary  importance  is  the  contouring  of  the 
airfoil  aft  of  the  leading  edge.  Thicker  contouring 
forward  of  "he  location  of  maximum  thickness  contributes 
to  delaying  flow  separation. 

3.  The  effect  of  increasing  the  pitching  rate  of  the  airfoil 
is  to  enhance  unsteady  lift  by  delaying  vortex  formation 
to  a  higher  angle  of  attack.  The  flow  solutions  present¬ 
ed  in  this  study  on  the  effects  of  reduced  frequency  are 
in  good  agreement  with  trends  observed  in  experimental 
work . 

4.  The  streamwise  pressure  gradient  required  for  flow 
separation  is  a  function  of  airfoil  leading  edge  radius 
and  is  independent  of  reduced  frequency  or  freestream 
speed . 


B.  RECOMMENDATIONS 

1.  Conduct  further  experimental  and  computational  investiga¬ 
tions  to  compare  Mach  number  effects  on  flow  development 
and  dynamic  stall.  The  investigations  should  be  conduct¬ 
ed  at  comparable  Mach  and  Reynolds  numbers. 

2.  Investigate  the  effect  of  airfoil  thickness  ratio  on 
dynamic  stall. 
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3.  The  understanding  of  harnessing  the  unsteady  lift 
generated  by  rapidly  pitching  airfoils  is  only  part  of 
the  knowledge  required  for  practical  implementation. 
Further  investigations  must  include  examining  pitching 
moment  development  and  how  to  minimize  the  adverse 
pitching  moments  produced. 

4.  Conduct  further  computational  work  at  higher  reduced 
frequencies  and  Mach  number  to  check  the  validity  of 
using  this  eddy-viscosity  model  at  those  flow  conditions. 
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APPENDIX  A  -  USING  THE  PROGRAM 


A.  THE  NAV I ER- STOKES  CODE  FOR  GENERATING  FLOW  SOLUTIONS 

The  main  program  NSE2D  (listed  in  Appendix  B  for  refer¬ 
ence)  reads  input  data  as  supplied  by  the  input  file  RSE.ld, 
the  grid  files  fort. 11  or  fort. 21  as  appropriate,  and  the  flow 
file  fort. 31.  The  program  forms  the  metrics  and  the  Jacobian 
by  calling  the  METRIC  subroutine.  After  calling  the  SLPS 
subroutine  to  perform  the  computations,  the  main  program 
outputs  the  rotated  grid  and  the  flow  solution  to  the  file 
fn  t.  20,  and  outputs  the  data  of  lift,  drag,  moment  and 
pressure  coefficients,  time,  angle  of  attack,  and  rotatjunal 
frequency  to  the  loads  file  fort. 3. 

The  subroutine  SLPS  is  the  primary  subroutine,  and  it 
calls  other  subroutines  to  perform  the  steps  in  the  ADI 
algorithm.  The  major  subroutines  that  SLPS  calls  are: 

Subront  ine _ Function _ 

MATRIX 1 , MATRIX2  Form  block  tridiagonal  matrices  for  the  £ 

and  T)  directions 

AMAT1,  AMAT2  compute  coefficient  matrices  c9F/(9£  and 

c5G/ 


STRESS 

DISSIP 


EXPBC 


supplies  the  viscous  terms 

adds  fourth-order  dissipation  terms  to 
the  right-hand  side  of  the  Navier-Stokes 
Equations 

applies  boundary  conditions 


113 


RES1 


computes  residuals  from  inviscid  part  of 
equations 

EDDY  applies  Baldwin-Lomax  mode] 

B.  DIRECTIONS  FOR  EDITING  INPUT  FILES  FOR  THE  PURPOSE  OF 
GENERATING  FLOW  SOLUTIONS 
1.  Generating  a  Grid 

Initially,  a  grid  must  be  generated  as  described  in 
Chapter  III  in  order  to  compute  the  flow  about  the  airfoil  in 
question.  For  this,  it  is  necessary  to  use  the  AIRFGR.F 
program.  With  the  rectangular  coordinates  of  the  desired 
airfoil,  edit  the  input  file  AIRFGR.IN  using  the  "ex"  or  the 
"vi"  editors.  Ensure  that  enough  points  are  defined  in 
regions  of  high  curvature  such  as  the  leading  edge  regions. 
Then  ensure  that  no  fort.  21  file  exists  by  renaming  or 
deleting  it.  At  this  time,  run  the  grid  program  by  typing 
AIRFGR  <  AIRFGR.IN 

After  the  program  has  run,  ensure  that  the  grid  is  satisfacto¬ 
ry  by  running  FL0T3D.  It  may  be  necessary  to  define  more 
coordinate  points  to  ensure  leading  edge  radius  definition, 
upper  or  lower  surface  curvature,  or  trailing  edge  closure. 
Accuracy  of  trailing  edge  closure  will  determine  wake  thick¬ 
ness.  The  AIRFGR.F  program  is  listed  in  Appendix  C  for 
reference . 

?.  Rotating  a  Grid 

The  grid  generated  in  part  1  is  for  a  zero  degree  angle 
of  attack.  Often  it  is  desired  to  compute  a  steady  flow  or  to 
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start  an  unsteady  flow  at  an  AOA  other  than  zero.  In  order  to 
do  this  the  ROTGR.F  program  must  be  used.  It  is  necessary  to 
edit  the  ROTGR.F  program  and  set  the  desired  rotated  grid 
angle  by  changing  the  ol  value  in  line  10.  After  the  program 
hap  been  edited,  then: 

1.  Compile  ROTGR.F’  by  typing  "cf77  rotgr.f". 

7.  Compiling  ROTGR.F  results  in  an  output  file  named  a. out. 
Move  a. out  to  rotgr  by  typing  "mv  a. out  rotgr". 

3.  Execute  the  program  by  typing  "rotgr". 

The  output  file  produced  by  rotgr  is  fort.  12.  It 
contains  the  desired  grid  rotation.  Copy  fort. IP  to 
fort. 11  or  fort. 21  as  necessary  for  future  flow  solu- 
t i one . 

The  commands  within  the  quotation  marks  should  be  typed,  not 
t';o  quotation  marks  themselves.  Upper  and  lower  case  are  not 
important . 

3.  Steady-State  Solutions 

Ensure  that  the  grid  is  rotated  to  the  desired  angle 
of  attack  of  the  flow  solution  (part  2).  The  rotated  grid 
.should  be  copied  to  the  file  fort.  21.  The  fort.  21  file  must 
have  the  grid  rotated  to  the  desired  angle  of  attack  tor  the 
steady-state  solution! 

Modify  the  NSE.IN  file  as  follows: 

line  2  : 

-  set  desired  computational  time  interval  ( DT ) 

-  set  ALFA,  ALFA I ,  ALFA1  to  zero 

-  set  REDFRE  to  desired  reduced  frequency 

-  set  AMINF  to  desired  Mach  number 

line  8  : 

-  set  number  of  time  steps  (ISTP);  usually  3000  are  neces¬ 
sary  for  convergence  at  A0A#0.  In, order  to  run  interac- 
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tively  and  not  exceed  time  limit,  set  ISTP=1000  and  run 
three  times. 

line  10: 

-  ensure  that  XREF=0.25 

-  set  TSHIFT=0 . 0 

-  set  REYREF  to  desired  Reynolds  number 
line  14: 

-  For  initial  run,  set  RESTART  OSCIL  RAMP  as  FALSE  TRUE 
TRUE.  To  restart  (2nd  and  3rd  runs),  set  TRUE  TRUE  FALSE 

line  2 2 : 

-  Set  time  to  0.0  for  1st  run.  Set  to  preceding  run  time  for 
restart . 


To  run  interactively  and  observe  output  as  it  develops,  type: 

NSE  <  NSE.IN 

To  run  in  the  background ,  type : 

NSE  <  NSE.IN  >  NSE. OUT 

When  the  run  is  finished  the  last  time  step  solution  will 
consist  of  grid  and  flow  solutions  in  file  fort. 20.  It  will 
be  necessary  to  separate  the  grid  and  flow  solution  to  restart 
the  second  run.  The  command  "SPLIT20"  will  move  the  grid 
solution  to  fort. 11  and  the  flow  solution  to  fort. 31,  and  will 
display  the  final  run  non-dimensional  time.  When  restarting 
the  second  run,  modify  NSE.IN  line  14  to  RESTART=TRUE .  A 
restarted  solution  uses  grid  file  fort.  11  and  flow  file 
fort .  31 .  Set  time  in  NSE.IN  to  time  specified  for  the 
previous  run  when  SPLIT20  was  executed.  Execute  the  next  run 
by  the  desired  method  as  before. 

After  convergence,  save  the  grid  and  flow  files 
( SPLIT20  to  fort. 11  and  fort. 31). 
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4.  Unsteady  Solutions 

Both  oscillatory  and  ramp  (constant  pitch  rate) 
solutions  can  be  obtained  by  the  NSE.F  program.  One  must 
ensure  that  the  desired  steady  state  initial  angle  of  attack 
grid  and  flow  solutions  are  obtained.  The  input  file  NSE.Ti.'  is 
modified  as  follows: 
line  2 : 

-  set  time  interval,  freestream  Mach  number  and  reduced 

frequency,  where  k=ic/2U(e  for  ramp  solution  and 

k=o>c/2Uai  for  oscillatory  solution. 

-  For  ramp  response:  set  ALFAI  to  steady-state  in.it  lal 
conditions;  ALFA,  ALFAI  are  unused. 

-•  For  oscillatory  response  of  form  a=a  «-*-a  .sin  ( <ot )  ,  set 

ALFA~a-,  ALFAl^a-,  ALFAI=  ( a  ^-a .)  since  response  will  start 

at  min  « . 

J.  i  ne  8  : 

-  set  I STP  6000-12000  depending  on  length  of  response 

desired . 

line  10: 

-  set  REYREF  to  desired  Reynolds  number. 

-  set  XREF=0.25 

net  TSHIFT  -  -0.5  so  that  response  starts  at  min  a  which 
is  1/4  cycle  ( -n/2  time  shift  as  o=a  ■,  +  a-  sin  ( ot~7t/ 2 )  ) 

from  mean  a. 


line  14: 

-  set  RE-START  to  TRUE 

-  set  OSCIL  to  TRUE  for  oscillatory  response  (else  FALSE) 

-  set  RAMP  to  TRUE  for  ramp  response  (else  FALSE) 

line  19: 

-  Program  flow  solution  outputs  can  be  selected  for  various 

angles  of  attack  desired.  Select  desired  A0A  solutions 
by  listing  them  multiplied  by  100  (i.e.  11  degrees  as 

1100).  The  outputs  are  listed  as  files  fort. 61  through 
fort. 72. 


line  22. 

-  Set  time  to  0.0  for  initial  run.  For  restart  set  time  to 
the  for t .  20  output  from  the  last  timestep  of  the  last  run. 
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Run  the  program  by  submitting  it  as  a  job  to  the  queue 
by  typing  "qsub  -It  7200  subunst".  The  subunst  file  is  a 
command  file  which  will  order  the  steps  for  the  Cray  when  the 
job  number  comes  up  in  the  queue.  With  at  least  1  1/2  hours 
in  computation,  go  run  5  miles,  eat  lunch,  or  play  with  the 
kids.  Upon  returning,  check  on  queue  status  by  typing  "mqs" . 
If  no  entries,  the  run  is  finished  and  you  can  split  the 
output  files  and  continue,  if  desired. 

Any  of  the  output  files  fort.  61  through  fort.  72  and 
fort. 20  can  be  split  into  grid  solutions  and  flow  solutions 
for  further  analysis.  Fort .  20  file  must  be  split  to  continue 
the  response  at  further  angles  of  attack,  and  to  use  the-  run 
time  as  the  next  run  NSE.IN  input  in  line  22.  For  continuous 
ramp  or  oscillatory  responses,  it  is  not  necessary  to  reset 
ALFA,  ALFAI ,  ALFA1 ,  or  RESTART;  it  is  only  necessary  to  reset 
desired  AOA  solutions  to  be  recorded  (line  19)  and  the  time. 

5 .  Example 

As  an  example,  to  run  an  oscillatory  solution  for 
a-1 Q±7sin ( t )  ,  it  is  necessary  to  get  a  steady  flow  solution  at 
a  =  3°  (3  =  10-7).  The  grid  developed  for  the  airfoil  must  be 
rotated  to  3°  using  the  ROTGR.F  program  and  input  to  fort  2.1. 
At  this  time  a  steady-state  solution  is  obtained  by  modifying 
NSE . IN  with  the  following  inputs: 
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ALFA 

ALFAI 

ALFA1 

REDFRE 

AMINF 

TIME 

RESTART  OSCIL  RAMP 

XREF 

TSHIFT 


0.0 
0.0 
0 . 0 

as  desired 
as  desired 
0 . 0 

FALSE  TRUE  FALSE 
0 . 25 
0 . 0 


Run  the  program  by  typing  "  NSE  <  NSE.IN  (  >NSE.0UT)  " 

After  1000  timesteps,  issue  the  command  SPLIT20,  so  that  the 
files  fort. 11  and  fort. 31  will  be  opened.  Edit  NSE.IN  with 
the  SPLIT20  time  and  set  RESTART=TRUE .  Run  the  program  again, 
and  then  a  third  time,  ensuring  convergence  (in  lift  coeffi¬ 
cient,  drag  coefficient,  etc.)  Save  the  final  fort. 11  and 
for t . 31  files. 

Then,  modify  NSE.IN  as  follows 


ALFA 

ALFAI 

ALFAI 

ISTP 

TIME 


10.0 
3 . 0 
7 . 0 

6000-12000 
0 . 0 


RESTART  OSCIL  RAMP  TRUE  TRUE  FALSE 

XREF  0.25 

TSHIFT  -0.5 


iai  .  ia2,  .  .  .  ial2  as  desired 

Run  the  program  by  issuing  the  command  "qsub  -It  7200 
subunst".  After  completion,  save  the  output  files  (fort.  61 
through  fort. 72)  and  fort. 20.  Do  a  SPLIT20 ,  note  time  and 
A0A ,  then  input  TIME  into  NSE.IN  for  further  runs. 

Graphs  of  the  output  files  may  be  obtained  by  using  the 
plot  program  PL0T3D.  Observe  output  files  (for  example 
fort.  65)  and  graph  by  typing  the  commands:  "split"  and  then  on 
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the  next  line  "65".  Outputs  will  be  in  files  fort. 11  (grid) 
and  fort. 31  (flow)  which  are  read  into  the  PLOT3L  input. 

A  typical  NSE.IN  file  for  the  unsteady  solution  of  the 
example  problem  at  k=0.02  and  M=0.4  would  be  as  follows. 


IMAX 

KMAX 

DT 

ww 

ALFA 

ALFA1 

ALFAI  REDFRE 

AMINF 

157 

58 

0.005 

2  . 

00  10.00 

7.00 

3.00  0.02 

0 .40 

ISPEC 

O 

(FLAG 

FOR  CHOOSING  DIFFERENT  SPECTRAL  RADIUS) 

3 

WW2X , 

WW2Y, 

WW4X,  WW4Y 

(EXPLICIT 

DISSIPATION  COEF.  FOR  X 

AND  Y) 

0.00 

0.00 

0.030 

0.030 

I STP 

NPER 

NOUT 

RES 

STRUNST 

9000. 

18000 . 

1000  . 

100  . 

0. 

REYREF 

DMIN 

XREF 

TSHIFT 

4.00 

0.00002 

0.25 

-0.5 

TSTAR1 

-1  . 

0 

RESTART, MULTIGRID  OPTIONS  SPECIFIED  IN  THE  NEXT  CARD 
TRUE  TRUE  FALSE 

CIRCOR  (  CIRCULATION  CORRECTION) 

TRUE 

31  127 

ITEL  I TEU 

0300  0400  0500  0600  0700  0800  0900  1000  1100  1200  1300  1400 
ial  ia2  ia3  ia4  ia5  ia6  ia7  ia8  ia9  ialO  iall  ial2 
TIME 
0.0 
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SOURCE  PROCRAM 


1 1:4&04  am 


SOURCE  TEXT 


ALMEAN  -  ALFA 

ALTAI  -  A1FA1  •PI  /  1B0. 

ALFAMAX  -  ALTAI  *►  ALMEAN 

ALFACR  •  ALFAMAX 

ALFAS  IS  THE  ANGLE  OF  AIRFOIL  HITS  RESPECT  TO  TREE STREAM  AT  STEADY-STATE 
OR  INITIAL  POSITION  OF  UNSTEADY  MOTION. 

IT  SHOULD  BE  SET  ACCORDING  TO  THE  TYPE  OF  MOTION 
ALFAS  -  ALMEAN  -  ALFAI *COS <  0 .  ) 

UINF  -  AMINF  •  COS (ALFA) 

VINF  -  AMINF  •  SIN ( ALFA) 
uuif  -  amint 
▼iof  -  0.0 

call  rotgxid<  laax,kjaax,  alf  a  ) 

DO  7  1-1, IMAX 
DO  7  K-l, KMAX 
Q1(I,I)*1. 

Q2(I,I)-UINF 
Q3(I ,K)»VINF 
Q4(l,K)-TOTEN 
7  CONTINUE 

IF(RSTST)  THEN 
REWIND  11 

READ  (11)  IMAX  ,  KMAX 

READ  (11)  (  (  X  ( 1  ,  K  )  ,  I-l.IMAX  ),  f-l.KMAX  ) , 

4  (  (  Z(I,K),  1-1, IMAX  ),  K-l ,  KMAX  ) 

rewind  31 

READ  (31)  IMAX  ,  IMAX 

READ  (31)  FSMACH ,  ALFAD,  REREAL,  TIKE 

READ  (31)  <  (  Q1(1,K),  1-1, IMAX  ),  X-1,KMAX  )  , 

4  (  (  Q2(I,X),  1-1, IMAX  ),  K-l , IMAX  ), 

4  (  4  03{I,K),  I- 1, IMAX  ),  K-l , IMAX  ), 

4  (  (  Q4 ( I , K ) ,  1-1, IMAX  ),  K-l , IMAX  ) 


if(  tiae  ,qt.  500.  )  ti»e  -  0. 
if(  alfad  . eq  0  )  ti»e  -  0. 

TS TART  -  TIME 

I F ( TSTART . GE . 0 . )  TIKE  -  TSTAFT 
IF( .NOT. (RSTRT))  TIME  -  0. 

NR  I TE  ( 6 , 6  8  )  TIME 

68  FORMAT (/, 4 0HTHE  CALCULATIONS  STARTED  AT  TIKE  T  •  ,F12.«,/) 

ASTART  -  ALFA  4  ALFAI  •  SIN(  2 ‘REDFRE 'AMINF  «  TIME  4  TSHIFT  ) 
ASTART  -  ASTART  •  <  180.  /  PI  ) 

I DONE  “  TIME  /  DT 

WRITE( 6,681)  ASTART,  I DONE 

681  FORMAT^ 2X . 1 0HALFA5TART* ,  F2  0 . 4 , 9X ,  2 5HX TERATIONS  DONE, 5X, 16,/) 
CALL  METRIC 

IKIN  “  (  ITEU  -  I  TEL  )  /  2  ♦  I  TEL 

I  LON  -  2  *  I  MIN 

CHD  -  X( I  TEL, 1 )  -  X ( ININ , 1 ) 

NSTPT  -  NSTP  4  NSTPP 


->  STARTING  FROM  A  STEADY- SATE  SOLUTION  AT  CERTAIN  ANGLE  <Ao> 
GIVE  A  TIME  SHIFT  <TSHIFT>  TO  START  FROM  AO  AND  SET  INITIALY 
TIME-0.,  KEEP  TSHIFT  THE  SAME  THROUGH  THE  UNST .  CALCUL. 


IF  (  PITCH  )  THEN 
DTP  -  PI  /  (NPER* REDFRE* AMINF) 

dt  •  dtp 

end  if 

alfacr  -  pi  /  9. 
alfamax--  alaean  4  alfal 
IF C  RAMP  )  READ  (5,*)  TIME 
DO  1000  JTN*1 ,NSTP 
TIME  -  TIME  ^  DT 


REDUCE  THE  TIME  STEP  TO  HALF  AT  THE  PEAK  OF  THE  CYCLE 


IT  (PITCH)  THEN 


IF(  ALFA  GT.  ALFACR  )  THEN 

DT  -  DTP  •  |  1.  -  . 5* (ALFA- ALFACR )/( ALFAMAX- ALFACR)  ) 

ELSE 
DT  -  DTP 
END  IF 

OMEGA  «  2.  *REDFRE*AMINF  •  COS(  REDFRE*  2  . ‘TIME* AMINF  TSHIFT) 
1  * ALFAI 

ALOLD  »  ALMEAN  +  ALFAI  •  SIN<  2.  •  REDFRE  *  AMINF  * 

1 (TIME  -  DT)  4  TSHIFT  ) 

ALFA  •  ALMEAN  4  ALFAI  •  SIN(  REDFRF  *  2.  *  TIKE  •  AMINF 
1  4  TSHIFT  ) 

ALFAD  *  ALFA  *45/  ATAN(1.) 

DALFA  •  ALFA  -  ALOLD 

CALL  ROTGRID(  IMAX.  KMAX,  DALFA) 

CALL  METRIC 
END  ir 

IF  (RAMP)  THEN 

OMEGA  -  2  *  REDFRE  •  AMINF 

ALOLD  -  OMEGA  *  (  TIME  -  DT  ) 

ALFA  •  OMEGA  •  TIME 

ALFAD  -  Air A  *45./  ATAN ( 1 . ) 

DALFA  -  ALFA  -  ALOLD 

CALL  ROTGRID( IMAX, KMAX, DALFA) 

CALL  METRIC 

end  xr 

ALFAD  -  ALFA  *45/  ATAN ( 1 . ) 


CALL  SLPS(ITN.ISPEC) 

APPLICATION  Of  BOUNDARY  CONDITIONS 


CALL  EXPBC(CL) 

CALL  LOAD ( CL, CDP, CDF, CM, ALFAS, XFEF) 


PRINT  OUT  PRESSURE  AT  THE  SURFACE 


IF( I  IN/ 50 *50 . EQ. ITN )  THEN 
WRITE  (6,19) 

WRITE  (6,13)  ITN,  TIME  ,  DT 

IF( PITCH  OR.  RAMP)  WRITE  (6,3500)  ALFAD, OMEGA 

lF{ITN/50''  *  500.  EQ.  ITN)  CALL  CPPLOT(  ITEL,  ITEU  ,  AMINF  ,X  ( 1 , 1 ),  Z(1 , 1 )  } 
IF(ITN/40'J*400.EQ  ITN)  WRITE(i,12)  ( I  ,CF  ( I  )  ,1-1 ,  IMAX  ) 

WRITE  ( 6 , j 000 )  CL  ,  CDP  ,  CDF  ,  CM 
END  IF 

TIKPI  -  TIME  /  AMINE 

JF( ITN/50* 50  EO  ITN) WRITE (3 , 8000 ) TIKE , ALFAD , OMEGA , CL, CDP , CDF , CM 
IF(  I TO/1 000*1 000. EQ. I TO)  THEN 
DO  100  T  •  IKIN  ,  ITEU 


SOURCE  PROCRAM 

DATE  1/20/90 

nse2d.f 

T,ME  1 1:46.-04  am 

SOURCE  TEXT 


C  XOC  -  <  X(I,1)  ~  J  )  /  CBD 

C  niTC(3/7000)  XOC  ,  PSOR(  I LOW— I)  ,  PSOR(I) 

C  WTTE(  9, 7000 )  XOC  ,  CF(ILOW-I)  /  CF(I) 

C  100  CCWTIWUE 

c  nrotr 

CWO  TO  CALCULATE  THE  AVERAGE  VALDES  Of  CL,  CD,  CM 
5CL  -  SCL  ♦  CL 
SCDF  -  SCDF  ■»  CDF 
SCDP  -  SCDP  +  COP 
SCM  «  SCM  •*  CM 

C«r  TO  NEXT  cw  for  user  specified  plots 

C  DO  $84  1  -  1  ,  IMAX 

C  DO  984  I  -  1  ,  I MAX 

C  *84  D*10< I , X )  -  ABS {  QL(1,R)  -  DRflOJI,!)  ) 

IRES  -  IFIX(RES) 

IF(ITN/IRES*IRES.EQ.ITN)  THEN 
IPLOT  -  (  NSTPP  +  IW  )  /  IRES 
R£SD( IPLOT)  -  0. 

CLH{ IPLCT )  «  CL 

CDPH( IPLOT)  •  CDP 
DO  985  R  «  1  ,  KMAX 
DO  985  1*1,  TMAX 

C  985  RESD  ( IP  LOT )  *  AKAXl(RESDOPLOT)  #DRHO(  I ,  K  )  > 
985  RESD( IPLOT )  -  RESDL2 ( ITN } 

ENDIF 


IF(  ITN/NOt'T*NOUT  .EQ-  ITN  )  THEN 
ir (  ITN  .EQ.  1 *NOCT  )  IOUT  -  31 
IF<  ITN  .FQ.  2*NOUT  )  IOUT  -  32 


IF(  ITN  . EQ. 
IF<  ITN  .EQ. 


5*NOUT  )  IOUT 
6*NOUT  )  IOUT 


IF(  ITN  .EQ.  2*NOUT  )  IOUT  -  32 

ir<  ITN  . EQ.  3-NODT  )  IOUT  -  33 

I F (  ITN  .EQ  4  *NOUT  )  IOUT  -  34 

IF(  ITN  .EQ.  5*NOUT  )  IOUT  «  35 

IF<  ITN  .EQ.  6*NOUT  )  IOUT  -  36 

IF(  ITN  . FQ.  7 *NOUT  )  IOUT  -  37 

IF(  ITN  .FQ.  8*NOUT  )  IOUT  -  38 

IF(  ITN  .EQ.  9-NOUT  )  IOUT  -  39 

IF(  ITN  .EQ.  10*NOUT  )  IOUT  -  40 

IF(  ITN  -EQ.  ll*NOUT  J  IOUT  -  41 

IT(  ITN  ,EQ  12«NOUT  )  IOUT  «  42 

IF<  ITN  .EQ.  1 3 *NOUT  )  IOUT  *  43 

IF{  ITN  .EQ.  14  *NOUT  )  IOUT  -  44 

IF (  ITN  .EQ.  15*NOUT  )  IOUT  «  45 

IF<  ITN  -EQ.  16-NOUT  )  IOUT  -  46 


iout  *  20 
REWIND  IOUT 

WRITE  (IOUT)  I MAX  ,  KMAX 

WRITE  (IOUT)  (  (  X ( I , K ) ,  1*1,1 MAX  ),  K-l , KMAX  ) 

4  (  (  Z( I , K  )  ,  I-1,IMAX  ),  K-l, KMAX  ) 

WRITE  (IOUT)  IMAX  ,  KMAX 
WRITE  (IOUT)  AMINF ,  ALFAD,  REREAL,  TIKE 
WRITE  (IOUT)  (  (  Ql ( 1 , K ) ,  1-1, IMAX  ),  K«1,KMAX  ) 

i  (<  02(1,1),  1-1,1  MAX  ).  K  =»  1  ,  KMAX  ), 

4  (  (  Q3(I,K),  1-1 , IMAX  ),  K-l.KMAX  ) 

(  (  04 ( I , K ) ,  1-1 , IMAX  ),  K«1 , KMAX  ) 


REWIND  IOUT 


IF( 

I ALFAD 

•  EQ. 

lal 

.OF. 

I ALFAD 

•  EQ. 

U2 

OR 

3 

I ALFAD 

■  EQ. 

i«3 

.OR. 

3 

I ALFAD 

•EQ. 

1*4 

.OF. 

% 

IALFAL 

EO- 

1*5 

.OR. 

3 

IALFAD 

•  EQ. 

1*6 

.OF. 

3 

IALFAD 

.EQ. 

1*7 

.OR. 

3 

IALFAD 

•  EO 

1*8 

.OR. 

3 

IALFAD 

.EQ. 

1*9 

.OF. 

3 

IALFAD 

-EQ. 

1*10 

OR. 

3 

IALFAD 

.EQ. 

1*11 

.OR. 

3 

IALFAD 

■  EQ. 

1*12 

)  THEN 

XF( 

IALFAD 

•  EQ. 

1*1 

)  IAOUT 

61 

IF( 

IALFAD 

.EQ. 

1*2 

)  IAOUT 

62 

IF( 

IALFAD 

.EQ. 

1*3 

)  IAOUT 

63 

XF( 

IALFAD 

•  EO. 

1*4 

)  IAOUT 

64 

IF< 

IALFAD 

.EO. 

1*5 

)  IAOUT 

65 

1F( 

IALFAD 

.EQ 

1*6 

)  IAOUT 

66 

IF( 

IALFAD 

•  EQ. 

1*7 

)  IAOUT 

67 

XF( 

IALFAD 

.eq. 

1*6 

)  IAOUT 

68 

IF( 

IALFAD 

.EQ. 

1*9 

)  IAOUT 

69 

IF( 

IALFAD 

.EQ. 

1*10 

)  IAOUT 

70 

ir( 

IALFAD 

.EO. 

1*11 

)  IAOUT 

71 

XF( 

IALFAD 

eq. 

1*12 

)  IAOUT 

72 

REWIND  IAOUT 

WRITE  (IAOUT)  IMAX  ,  IMAX 

WRITE  (IAOUT)  (  (  1(1,1),  I-l.IKAX  ),  1-1, KMAX  ), 

4  (  (  2(1.1),  I - 1 , IMAX  ),  1-1, IMAX  ) 

WRITE  (IAOUT)  IMAX  ,  KMAX 
WRITE  (IAOUT)  AMINF,  ALFAD,  REREAL,  TIME 
WRITE  (IAOUT)  (  (  01(1. I).  I- 1 , IMAX  ),  1-1, KMAX  ) 

4  (  (  02(1,1),  I - 1 , IMAX  ),  1-1, KMAX  ), 

K* 1 , KMAX 


REWIND  IACUT 


(  (  03(I,D  , 
(  (  04(1,1), 


1-1 , IMAX  ) , 
1*1, IMAX  ), 


1000  CONTINUE 
C 

REWIND  20 

WRITE  (20)  IMAX  ,  KMAX 

WRITE  (20)  <  (  X(I.K),  1*1, IMAX  ),  1-1, IMAX 

4  (  (  2(1, K),  1*1,1 MAX  ),  K* 1 , KMAX 

WRITE  (20)  IMAX  ,  KMAX 

WRITE  (20)  AMINF,  ALFAD.  REREAL,  TIME 
WRITE  (20)  (  (  Ol(l.K),  1*1, IMAX  ),  K“  1 , KMAX 
4  (  (  02(1.1),  I -1, IMAX  ),  R* 1 , IMAX 

4  (  (  03(1, K),  1-1 , IMAX  ),  K* 1 , KMAX 

4  ((  04(1.1),  1*1.1 MAX  ),  K= 1 , KMAX 


SOURCE  PROCRAM 

nse2d.f 


1/20/90  PACt* 


T,ME  11^04am 


SOURCE  TEXT 


C 

alfat  -  olta  •  (180/pi) 
vrite(6,667)  alfaf,  ti*e 

467  torwat ( 17hANGLE  OF  ATTACK  - ,F15 . 10, 5X, 5HTIKE- , F10 . 4 ,//) 

C 

c  primt  out  velocity  motile 

c 

WRITE( 6,668) 

668  FORMAT (// , 14X , 2HUX ,141, 2HUZ ,  10X , 8BDB4SITY  , 6X , 6HPRESSUFE ,// ) 

CMIIL  •  AMXNF  /  RFYRXF 
DO  4000  I  -  70,130,2 
S  -  0. 

DO  4000  F  -  2  ,20 

S  -  S  ♦  S0RT(  (X(I,JC)-X(I,E-1))*«2  +  ) 

V  -  (  Q2 ( I , K )/Ql ( 1 , K )  ) 

V  -  (  Q3 ( I , K }/Ql ( I , I )  ) 

UTOT  •  0*U  +  V*V 

P  -  (GAKMA-l.)M  04(1,1)  -  .  5*01  (I,I)*BTOT  ) 

WRITE  (6,2002)  I,K,  U,V  ,  Q1(I,R)  ,  P 
4000  CONTINUE 
C  CALL  OUTKVC  ( RERFAL ) 

STOP 

C  1  PORMAT( 12A6  J 

1  FORMAT(12A8) 

2  FORMAT (225, 7 FI 0.0) 

11  FORMAT (7F1 0.0) 

3  FORMAT (/, 5X , 5BIMAX* , I 1 0 ,// , 5X , 5HKMAX* ,110, 

1/,SX  58  DT-,F20.8,/,5X,5BWW  * , F20 . 8 , / , 51 , 5HWWI  -,F20.8, 

2/ , 5X , 5BALFA= , F20 . • ,/, 5X , 6BALFA1* , F20 . 8 , / , 5X , 6BALFAI* , F20 . 8 ,/, 

+  5X , 7HREDFRE* , F18 . 8 , / , 5X , 6HAM3NF- ,F2 0 . 8 , /, 5X , 5BXREF- , F21 . 8 ) 

4  FORMAT ( 1H1 . 5X , 10A6 ) 

12  FORMAT (8(14, F10. 4) ) 

19  FORMAT  ( /  ) 

C  22  FORMAT ( 2F1 0.6,15) 

3  3  FORMAT ( 5> , 5HI STP- , I  5 , 5X , 5HTIME- ,  F9 . 5 , 5X , 3HDT* , F9 . 5 ) 

2002  FORMAT (215, 5E 14. 6) 

3000  FORMAT ( 5X , 38C L=  ,  FI  0 . 4 , 5X , 4BCDP- , F10 . 4 , 5X , 4BCDF* , F10 .  4 , 5X  , 3HCM= , 
•*ri0.4) 

3  500  FORMAT ( 5X , 5BALFA- , F10 . 4 , SX , 6HOMEGA- , F10 . 4 , 5X , 2BH* , FI 0 . 4 , 5X , 5HHDOT* 
1  ,  F10 . 4  ) 

3700  FORMAT ( 5X , 7HREYREF- , F2 0.4) 

C5000  FORMAT( 21*10-6  ) 

6000  FORMAT (811 4- 8) 

7000  FORMAT  (31  16.**) 

8000  FORMAT ( 91 13.6) 

8  001  FORMAT  (3.'.  ,  8BCL(  AVG  )  «=  ,  FI  2 . 7 , 4X  ,  9HCDF ( AVG  )» , F12 . 7 , 4X ,  9HCDP  ( AVG  )s  , 

+  F12 . 7 , 4X  , 8BCM( AVG) * , F12 . 7) 


SOURCE  PROGRAM 

nse2d.f 


DATE 


TIME 


LINE  i 


SOURCE  TF^T 


SUBROUTINE  AMAT1 ( l , IKX1 , XIX , XIZ.XIT) 


SUBROUTINE  AMAT1 


PARAMETER  ( IX-1 80 , fcx-60 ) 

C0(«0N/rU)N/01(IX,RX)  ,Q2 ( IX , XX  )  , Q3 ( IX  , XX ) , Q4 { IX , XX ) 
C0MM0N/PERTR/D01(IX,XX) ,DQ2(IX(XX) ,DQ3 (IX, XX } , DQ4 ( IX , XX ) 
COHHON/AM  A( 4 , 4 , IX ) 

COMMON  /  p  AR/GAMKA ,  REYREF ,  ALFA ,  ALFA1 ,  REDFRE  ,  AMIN  F  ,  ALFAI 
DIMENSION  1IT( IX , XX ) , XIX ( IX , XX ) ,  XI  Z  ( I X  ,  XX  ) 

REAL  XG,X1,X2 

AMAT1  COMPUTES  THE  COEFFICIENT  MATRIX  DE/DQ  DURING  XI  SWEEP 


GM1 

DO  1000  I 

X0 
XI 
X  2 
V 
H 

EBYR 
PHI  2 
THETA 

A( 1 , 1 , 1 ) 
A(  1 «  2 . 1 ) 

w.y.i) 

*(1.4,1) 
*(2,1,1  ) 
A( 2 , 2 , 1 ) 
*(2,3,1) 
*(2,4,1) 

*(3,1,1) 

*(3,2,1) 

A( 3 , 3 , 1  ) 

*(3,4,1) 

*(4,1,1) 

*(4,2,1) 

*(4,3,1) 

*(4.4,1) 

CONTINUE 

RETURN 

END 


GAMMA  -  1. 

2  ,  IMX1 
XIT(J ,X) 

XIX(I ,K) 

XIZ(I,X) 

02(1, X)  /  Q1 ( I , K ) 

Q3(I,X)  /  Q1 ( I , X ) 

04 { i , r  >  /  qi { i . x ) 

0.5  •  CK1  *  (U  *  V  *  w  •  K) 

XI  *  U  +  X2  •  V 
KO 
XI 
K2 
0 

XI  *  PHI  2  -  U  *  THETA 
HO  ■*  THETA  -  XI  •  (GM1  -  1.)  •  U 
X2  •  V  -  GM1  •  XI  *  N 
XI  •  GM1 

X2  *  PHI  2  -  N  *  THETA 
XI  *  N  -  X2  •  GM1  •  U 
FO  THETA  -  X2  *  (GM1  -  1.)  •  N 
X2  •  GM1 

THETA  •  (2.  *  PHI 2  -  GAMMA  •  EBYR) 

XI  •  (GAMMA  *  EBYR  -  PHI 2)  -  GM1  •  U  •  THETA 

K2  *  (GAMM.A  •  EBYR  -  PHI2)  -  GM1  *  V  •  THETA 

XO  •»  GAMMA  «  THETA 


1/20/90 

PACE# 

11:4004  am 

5 

SOURCE  PROCRAM 

nse2d.f 


date  1/20/90  f*CE  * 

™«  11:4&04  am 


SOURCE  TEXT 


SOURCE  PROCRAM 

DATE  1/20/90 

nse2d.f 

™E  11S4&04  am 

LINE  # 

499 

“SOO 

C« 

SOI 

c* 

S02 

C* 

S03 

c* 

504 

IaC 

c- 

SOURCE  TEXT 


SUBROUTINE  SLPS ( ITN , ISPEC) 


SUBROUTINE  SLPS 


PARAMETER  ( IX-180 , kx-60 ) 

CONNON/rLOK/01 (IX , XX ), Q2 < IX , XX ) , Q3 ( IX , XX ) , 04 <ZX , IX) 

COMMON'/ ri  X /OMEGA ,  HOOT 

COMMON/PER TR/DQ1  <  IX ,  KX  )  , DQ2  { IX ,  KX )  ,  DQ3  ( IX ,  KX )  ,  DQ4  ( IX  ,  KX ) 

COMMON/AH  M«  <4,  IX) 

COKMON/TR I Dl/DD ( 4 , 4 , IX , KX ) 

COMMON/TF: D2/MM ( 4 , 4 , IX , KX ) 

COMMON/TRI D2/EE( 4,4, IX  # XX) 

COMMON/TR I D4 /GG < 4 , I X , XX J 

COMMON /PAR /GAMMA , REYREF , ALFA , ALFAI , REDFRE , AMINF , ALFAI 
COMMON/DGR I D/ DT , I MAX , UMAX , I TE L , I TEU 
COMMON/GR I D/Y ACOB { I X , XX ) 

COMMON/D AVP/KK , WW I ,  NW2 X , WN2Y , WW4 X f NW4 Y 

COMMON/WTKIX/  XIX( IX , XX ) ,  II  Z.  { IX  ,  KX  ) , ZETAX ( IX , XX ) , ZETA2 ( IX , KX ) 

1  ,XIT( IX , KX ) , ZETAT< IX , KX ) 

COMMON/RADI  US /  SPECX ( IX , KX) , SPECY ( IX/ XX ) 

COMMON/ L2  NORM/  RESDL2 ( 10000 ) 

REAL  MM , DL LTAT ( Z  X , KX ) 

Tffi  SITBROL TINE  SLPS  DOES  TRE  BULK  OF  THE  WORE  FOR  THE  ADI  ALGORITHM. 
IT  CALLS  FLUX  AND  COMPUTES  TIGHT  HAND  SIDE  DURING  THE  TWO  SWEEPS, 
ASSEMBLES  THE  COEFFICIENT  HARICES,  ADDS  IMPLICIT  AND  EXPLICIT 
DISSIPATION  AND  CALLS  THE  TRIDIAGONAL  SOLVER  TO  OBTAIN  THE  FINAL 
SOLUTION. 

I*U  -  I  MAX  -  1 

I M2  •  I MAX  *  2 

KK1  «  KMAX  -  1 

EM 2  -  KMAX  -  2 


THE  DISSIPATION  TERMS  ARE  CONSTRUCTED  AND  STORED  IN  THE  ARRAYS  DO I , 
DQ2 , DQ3  AND  DQ4 . 

CALL  SPTCT ( ISPEC) 

CALL  DISS3P 

ON  TO  DQ1.DQ2,DQ3  AND  DQ4 . 

THE  RIGHT  HAND  SIDE  AT  KNOWN  TIME  LEVEL  IS  NOW  COMPUTED  AND  ADDED 
CALL  RESI 

DO  8999  X  -  2  ,  JW1 
lx,  8999  I  «  2  ,  IMI 
DELTATU  ,  K )  «=  '  . 

I F( ( REDFRL . LT . 0 . 001  )  . AND . ITS . EQ. 1 )  THEN 
DO  1  K  *  2  ,  KM1 
DO  11*2,  IMI 

DELTAT ( l , K )  -  1.0  /  (  X.  4  SQRT( ABS { YACOB< I , K ) )  )  ) 

ENDIF 


IF  VISCOUS  FLOW  IS  COMPUTED  THE  VISCOUS  TERMS  ARE  ADDED  TO  DQ1  ETC.  HERE 


!F(REYREF.'  0.)  CALL  STRESS(ITN) 
I— SWEEP . 

DTB  «  DT  •  0 . 5 

DTK  -  DT  •  WXI 

DO  3  T  *  2  ,  KM1 

CALL  AMATa { K , I MAX- 1 .XIX,XIZ,XIT) 


DO  4  II 
DO  4  12 
DO  5  I 

EE( I 1 , 12 , 1-1 , K ) 
DD( XI, 12, 1-1, K) 
CONTINUE 
CONTINUE 


,  4 
,  4 

,  IMI 

A(I1  ,12 ,1*1) 
A(I1, 12,1-1) 


IMPLICIT  DAMPING  ADDED  HERE. 


DO  6  I 

DTI 

DTWI 

DD(1 ,1 ,1-1 , K ) 
DD(2,2,1-1 , K ) 
DD  <  3 , 3 , 1 - 1 , K ) 
DD(4,4 ,1-1 ,K) 
EE(1 ,1 ,1-1  ,K) 
EE(2, 2,1-1 ,K) 
EE(3,3,I-1,X) 
EE ( 4 , 4 , I- 1 , K ) 
MK< 1 , 1 , 1-1 ,K) 
»«(2, 2,1-1 ,K) 
MH<3, 3,1-1 ,K) 
MM(4.4 ,1-1 ,K) 
CONTINUE 
CONTINUE 
DO  990  K 
DO  990  I 
GG( 1 ,1-1 , K  ) 
GG( 2 , I -1 , F ) 
GC( 3 , 1-1 , R ) 
GC(4,I-1,K) 
CONTINUE 


-  2  ,  IMI 

»  SPECX (I ,  X )  •  DTW  ' 
*  2.  •  SPECX ( I, K)  • 
DD ( 1 , 1 , I - 1 , K )  -  DTI 
DD(2,2 .1-1 ,K)  -  DTI 
DD( 3 , 3 , 1-1 , K )  -  DTI 
DD ( 4 , 4 , I - 1 , K )  -  DTI 
El.(1  ,1  ,1-1, K)  -  DTI 
EE ( 2 , 2 , Z-l , X )  -  DTI 
EE ( 3 , 3 , 1-1 , X )  -  DTI 
EE ( 4 , 4 ,1-1, K)  -  DTI 
1.  4  DTWI 
1.  +  DTWI 
1.  +  DTWI 
1.  4  DTWI 


■  2  ,  KM1 
»  2  ,  IMI 

-  DQ1 ( I , K )  •  DELTAT(I ,K) 

-  D02(I ,K)  •  DELTAT(1 ,K) 

-  DQ3(I,K)  *  DELTAT ( 1,1) 
«■  DQ4(I,XJ  *  DELTAT(I.K) 


DELTAT { I , K ) 
DELTAT ( I  ,  X) 


DELTAT( I , K ) 

YACOB ( I , K )  •  DTW* DELTAT { I ,  K ) 

*  Y  ACOB ( I  - 1 , K ) 

•  Y  ACOB ( I - 1 , K ) 

•  Y  ACOB ( I - 1 , K ) 

•  Y  ACOB ( I - 1 , K ) 

•  YACOB ( I 4l , K ) 

*  YACOB<l4l ,K) 

*  Y  ACOB ( I ♦ 1 , K ) 

*  YACOB(l4l,K) 


PERFORM  BLOCK  TRIDIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 

CALL  MATRX1 ( I MAX , KMAX ) 

DO  991  K  -  2  ,  KM1 

DO  991  I  -  2  ,  I Ml 

DQ1 ( I  , K )  •  GG(1  ,1-1, R) 

DQ2(I,K)  -GG(2,1-1,K) 

DQ3<I/X)  -  GG ( 3 , 1- 1 ,  K  ) 

D04  < I , K )  -  GG<4,1-1.K) 

CONTINUE 


K-SWEEP  BEGINS  HERE. 


DO  13  I  •  2  ,  IMI 

CALl  AMAT2( T , KMAX- 1 , ZETAX , ZETAZ , 2ETAT ) 


DO  15  II 
DO  15  12 
DO  15  K 

EE  ( 1 1 , 1 2  , 1  ,  K- 1 ) 

DD( II , 1 2 , 1 ,R- 1  ) 
CONTINUE 


I  ,  4 

1  ,  4 

2  ,  RM1 

A( 1 1 , 12 ,  K+l ) • DTH 
-A(I1,I2,K-1)»DTH 


DELTAT {1 ,K) 
DELTAT ( I , K ) 


STCONT  OR^ER  DAMPING  ADDED  HERE. 


SOURCE  PROCRAM 

nse2d.f 


PATE _ 1/20/90  PACE  * 

1U4&04  am  ® 


SOURCE  TEXT 


DO  16  K  -  2  ,  KN1 

DTI  -  SPECY(ZrK)  •  DTW  *  DELTAT(  I ,  K ) 

DTWI  -  2.  -  SPECY (1,1)  •  YACOB(I,K)  *  DTW*DELTAT( I ,  K) 

DD ( 1 , 1 , I , K-l )  -  DD( 1 , 1 , I , K-l )  -  DTI  •  YACOB( 1,1-1) 

DD(2,2,I,K-1)  -  DD ( 2 , 2 , 1 , K-l )  -  DTI  *  YAC0B(I,R-1) 

DD( 3 , 3 , I ,  K-l )  •  DD(3,3,X,K-Z)  -  DTI  *  YACOB(I,K-l) 

DD(4 ,4 , I ,K-1 )  ■  DD(  4 ,4,1, K-l)  -  DTI  •  YACOB(l,R-l) 

EE (1,1. I ,  K-l)  •  EE(1,1,I,K-1)  -  DTI  *  TACOB(I,R+l) 

EE( 2 , 2 , X ,  K-l )  »  EE ( 2 , 2 , I , K~1 )  -  DTI  •  YACOB ( I , *+ 1 ) 

EE( 3 , 3 , 1  ,  K-l )  -  EE(3,3 ,1 , K-l)  -  DTI  *  YACOB(I,K4l) 

EE( 4 , 4 , I ,K-1 )  -  EE<4 ,4,1 ,K-1)  -  DTI  *  YACOB(I,K+l) 

MM { 1 , 1 » I ,  K- 1 )  -  X.  4  DTWI 
MM ( 2 , 2 , X , K- 1 )  -  1.  +  DTWI 
MM(3,3,I,K-1)  -  1.  4  DTWI 
»<(4,4,I,K-1)  -  1.  4  DTWI 


EE( 2 , 2 , I  ,K-1 ) 
EE( 3 , 3 , 1 ,K-1 ) 
EE(4,4.X,R-1) 
Mt(l.l.Z,X-l) 

W(2,2,I,K-1) 

MM(3,3,I,K-1) 

16  M<(4,4,1,K-1> 
13  CONTINUE 

DO  17  K 
DO  17  I 
CG(1,I ,K-1) 

CG( 2 , I , K-l ) 

GG (3,1, K-l ) 
CG(  4 , X , K-l ) 

17  CONTINUE 


2  ,  KK1 
2  ,  I  Ml 
DQ1 <  I ,  K ) 
DQ2(I,K) 
DO 3 ( 1 ,  K) 
DQ4  (  I ,  K  ) 


PERFORM  BLOCK  TR I  DIAGONAL  MATRIX  INVERSION  FOR  THE  ENTIRE  PLANE 
CALL  KATRX2 ( I  MAX , KMAX ) 


DO  18  K 
DO  18  I 
DQ1(I,K) 
DQ2 ( I , K ) 
DQ3(I,K) 
DQ4 ( 1 , K ) 
18  CONTINUE 


2  ,  KK1 
2  ,  I  HI 
GG(l.Z.K-l) 
GG (2,1, K-l) 
GG( 3 ,1 , K-l ) 
GG (4*1, K-l) 


UPDATE  FLOW  VARIABLES  AT  INTERIOR  POINTS. 


967  CONTINUE 

RMAX 

0. 

RUMAX 

0. 

RVMAX 

0. 

EMAX 

0. 

DO  995  K 

2  ,  KK1 

DO  19  I 

2  ,  IM1 

01(1, K) 

* 

OKI  ,K) 

>♦ 

DQ1(I , K) 

*  YACOB ( I  ,  K ) 

02 ( I  ,  K  ) 

02(1, K) 

+ 

DO 2 ( 1 ,  K ) 

•  YACOB ( I , K ) 

03(1,1) 

03(1 , K ) 

DQ3 ( I , K) 

•  YACOB { I , K ) 

04(1, F) 

04(1, K) 

■* 

DQ4 ( I , K ) 

*  YACOB ( I , K ) 

19  CONTINUE 

DO  99b  I  *  2  ,  IM1 

IF  (RMAX . LT . ABS( DQ1 ( 1 , F) *YAC0B( I , K) ) )  THEN 
IR  *  I 

kr  *  r 

END  IF 

RKAX  «  AMAX1 (RMAX , ABS( DQ1 (I , K)  *  YACOB(I,K))) 

RUMAX  «  AMAX1  ( RUMAX  ,  ABS  ( DQ2  ( I  ,  K  )  *  YAC0B(1,K))) 

RVMAX  ■=  AKAX 1 ( RVMAX , ABS ( DQ3 ( I , K )  «  YACOB(l.K))) 

EMAX  -  AMAX1 ( EMAX , ABS ( DQ4 ( I , K )  «  YAC0B(I,K))) 

995  CONTINUE 


COMPUTE  L2  NORM  OF  Q1 ,  Q2 ,  Q3 ,  AND  Q4  RESIDUAL 


FUL2  -  0 
RVL2  »  0 
EL2  -  0 
DO  996  K  =  2  ,  KK1 

DO  996  I  *  2  ,  IK1 

RL2  •  RI2  ♦  D01 ( 1  .  F  )  •  *  2 
RUL2  -  RI  L2  *  DQ2  ( I  ,  K  )  •  •  2 
RVL2  -  RVL2  4  DQ3(1,K)*«2 
E 12  -  EL2  4  DQ4 ( I , F ) *  *2 

996  CONTINUE 

RL2  *  RL2  /  (  IM2*KK2  ) 

HU  12  -  Rl'L2  /  (  JM2  *  FK2  ) 

RVL2  «  PVL2  /  {  I M2  *  FM2  ) 

EL2  ■  EL2  /  (  I M2 • KM2  ) 

RL2  -  5C-HT(RL2) 

RUL2  »  SORT (RUL2 ) 

RVL2  -  SORT (RVL2 ) 

EL2  -  SORT (EL2 ) 

RESDL2(ITN)  -  RL2 

IF((ITN-1)/500»500.EQ  (ITN-1))  WRITE  (6,3002) 

IF(ITN/5P*50.EQ. ITN)  WRITE  (6,3001)  RMAX, RUMAX, RVMAX, EMAX , I R, KR 
IF(ITN/50*50.EQ. ITN)  WRITE  (6,6004)  RL2  ,FUL2  ,RVL2  ,EL2 
RETURN 

3002  FORMAT ( / , 5X , ’ DRMAX 1 , 1 1 X , ‘ DUMAX 1 ,11X, 1 DVMAX * , 11X , 1  DEMAX  * , 9X, 

1 'IR'  ,3X,  'KF' ) 

3001  FORMAT (4(E14.6,2X),2I5) 

6004  FORMAT* 4 (El  4. 8, 2X)  ) 

ENr 


SOURCE  PROCRAM 

DATE  1/20/90 

nse2d.f 

T,ME  11:4604  am 

SSEuii 


SOURCE  TEXT 


SUBROUTINE  MATRX1 (I MAX, KMAX) 


71.3 

C* 

SUBROUTINE  MATRX1 

* 

714 

c* 

* 

715 

c** 

• 

716 

PARAMETER  ( IX-180 , Ax-60 ) 

7)7 

COMMON/TRI Dl/DD ( 4 . 4 , IX , XX ) 

718 

COMHON/TRlD2/KM(4  ,  4  ,IX,KX) 

719 

COMMON/ TRI D3/EE ( 4 , 4 , IX , XX ) 

720 

COMMON/ TR I D4/  GG(4,IX,KX> 

72 1 

COMMON,  IZF.?.7/A('  ?,rX},»?<4 

4  77  ai  ,C(4 , 5, IX J 

722 

REAL  MM 

7» 

REAL  Lll ,  L2L  ,  L31 ,  L41 ,  L22  f  L32  ,  L42 , 133  ,  L4  3  ,  L44 

724 

2 , Lll ,  L2I  / L3I , L4I 

7|6 

c 

THIS  SUBROUTINE  PERFORMS  THE  BLOCK  TRI 01 AGONAL  MATRIX  I VERSION 

FOR 

7J7 

c 

AN  ENTIRE  PLANE  IMJRING  THE 

XI-  SHEEP 

72* 

729 

c 

KM1  -  KMAX  -  1 

750 

DO  I  11  -  1,  4 

751 

DO  1  K  -  2  /  RM1 

7S2 

AI  *  1 .  /  MM(1 ,1.1 ,K) 

735 

GG ( 11 , 1 , X )  -  GG( I 1,1,1)  *  AI 

734 

BB ( 1 1 , 1 , 1 ,  X )  •  EE(I1, l,l,K> 

•  AI 

735 

HB<I1,2,1,X)  -  EE { 11,2,1 , X  ) 

*  AI 

736 

HE ( 11 , 3 , 1 , K)  -  EE ( 11 , 3 , 1 ,K) 

*  AI 

7  37 

HB ( 1 1 , 4 , 1 , K)  -  EE (11, 4,1, K) 

•  AI 

738 

1  CONTINUE 

739 

740 

c 

DO  1000  I  -  2  ,  I MAX  -  2 

741 

DO  5  11-  1  ,  4 

742 

DO  2  X  -  2  ,  XM1 

743 

C(I1 ,1,K)  -  GG ( I 1 , I , K )  - 

DD( 1 1 , 1 , 1 , X )  *  GG( 1 , 1-1 ,  X  ) 

744 

1 

DD ( 11,2,1 , X )  •  GG(2,1-1,X) 

745 

2 

DD ( I 1 , 3 , I , X )  •  GG(3,I-1,X) 

746 

3 

DD ( I 1 , 4 , I ,X)  «  GG(4,I-1,X) 

747 

2  CONTINUE 

74* 

DO  5  12  -  1  ,  4 

749 

DO  5  X  -  2  ,  XK1 

750 

A(  1 1 , 1 2  ,  X )  -  MM ( I 1 ,12, I, X)  - 

DD ( 2 1 , 1  , 1 , X )  *  -HH(1, 12,1-1, X) 

751 

1 

DD ( I 1 , 2 , I , X )  •  HH ( 2 , I 2 ,1-1 , X ) 

752 

2 

DD(Il/3/I,K)  •  HH( 3 ,12 ,1-1 , X) 

753 

3 

DD ( I 1 , 4 , I , X )  *  HB ( 4 , 12 , 1-1 , X ) 

754 

C(I1,12*1,X)«  EE (11, 12, I, K) 

755 

5  CONTINUE 

756 

DO  3  X  -  2  ,  KM1 

757 

Lll  -  A( 1 , 1 , X ) 

758 

Lll  -  1.  /  Lll 

759 

U1  2  -  A( 1 , 2 , X)  •  Lll 

760 

U13  -  A ( 1 , 3  ,  X )  •  Lll 

761 

U1 4  -  A  (  1 , 4  ,  K  )  •  Lll 

762 

L21  -  A(  2 , 1 ,  X ) 

763 

L31  -  A( 3 , 1 , X ) 

764 

L41  -  A ( 4 , 1 , X ) 

765 

L22  -  A ( 2 , 2 , X  )  -  L21  *  U12 

766 

L2I  -  1.  /  L22 

767 

U23  -  ( A ( 2 , 3 , X )  -  L21  •  U13)  •  L2I 

768 

U24  -  ( A  ( 2  . 4  ,  X )  -  L21  •  U14)  -  L2I 

769 

L32  -  A(  3 , 2  ,  X  )  -  L31  *  U12 

770 

L42  -  A { 4 , 2 , X )  -  L41  •  U12 

771 

L33  -  A (  3 , 3  ,  X )  -  L31  *  U13 

-  L32  •  U23 

772 

L3I  -  1  .  /  L33 

773 

U34  -  ( A( 3 , 4 , X )  -  L31  •  U14 

-  L32  •  U24)  *  L3I 

774 

L43  -  A( 4 , 3 , X )  -  L41  •  U13 

-  L42  •  U23 

775 

L44  -  A ( 4 , 4 , X }  -  L41  •  U14 

-  L42  «  U24  -  L41  •  U34 

776 

L4I  -  1.  /  L44 

777 

C(1,1,X)  -  C(l.l.K)  *  Lll 

778 

C( 2 , 1  ,  X  )  -  (C(2,l ,X)  -  L21 

•  C  ( 1 , 1 ,  K )  )  •  L2I 

77  9 

C ( 3 , 1 ,X)  -  (C( 3 , 1 , X  )  -  L31 

•  C(1,1,X) 

780 

1  -  L32 

*  C  (  2  / 1  ,  X )  }  *  L3I 

78) 

C ( 4 , 1 , X  )  -  (C( 4 , 1 , X)  -  L41 

•  C ( 1 , 1 , X )  -  L42  •  C ( 2 , 1 , X ) 

782 

1 

-  i,4j  *  C  (  3 , 1 ,  X )  )  *  L4I 

783 

C( 1 , 2  ,  X  )  -  C( 1 , 2 , X  )  •  Lll 

784 

C ( 2 , 2 , X  )  -  (C( 2 , 2 , X)  -  L21 

*  C( 1 , 2 , X) )  «  L2I 

785 

C(  3 , 2  ,  X  )  -  (C (  3 , 2  ,  X)  -  L31 

•  C(1,2,K) 

786 

1  -  L32 

•  C ( 2 , 2 , X ) )  *  L3I 

787 

C ( 4 , 2  ,  X  )  •  (C( 4 , 2 , X)  -  L41 

•  C ( 1 , 2 , X )  -  L42  •  C ( 2 , 2 , X ) 

788 

1 

-  L4  3  •  C ( 3 , 2  ,  X )  }  •  L4I 

789 

C(1 ,3,X)  -  C(1 ,3,X>  *  Lll 

790 

C( 2 , 3  ,  X  >  -  (C( 2 , 3 , X)  -  L21 

•  C  ( 1 , 3  ,  X ) )  •  L2I 

79 1 

C(  3 , 3  ,  X)  -  (C(  3 , 3  ,  X)  -  L31 

-  C( 1 , 3 , X ) 

792 

1  -  L32 

•  C  (  2 , 3  ,  X  )  )  •  L31 

293 

C( 4 , 3 , X )  -  (C( 4 , 3 , X )  -  L41 

•  C ( 1 , 3 , X )  -  L42  •  C( 2 , 3 , X) 

794 

1 

-  L43  •  C  (  3 , 3  ,  X  )  )  •  L4I 

79S 

C( 1 , 4  ,  X  )  -  C( 1 , 4 , X )  •  Lll 

796 

C( 2 , 4 , X  )  -  (C( 2 ,4 , X)  -  L21 

*  C ( 1 , 4  ,  X } )  •  L2I 

797 

C( 3 , 4  ,  X  )  -  (C( 3 , 4 , X )  -  L31 

•  C  ( 1 , 4  ,  X  ) 

798 

1  -  L32 

*  C( 2 , 4  ,  X ) }  •  L3 1 

799 

C(4 ,4 , X )  -  (C( 4 , 4 , X )  -  L41 

•  C ( 1 , 4 , X }  -  L42  •  C ( 2 , 4 , X ) 

800 

1 

-  L4  3  •  C  (  3 , 4  ,  X )  )  •  L4I 

801 

C{ 1 / S  <  X)  -  C( 1 , 5 , X)  «  Lll 

*02 

C(2,5,X)  -  (C(2,5,X)  -  L21 

•  C(  1 , 5,  X  ) )  *  L2I 

803 

C( 3 , 5 , X )  -  (C(3,5,X)  -  L31 

•  C  ( 1 , 5  ,  X  ) 

804 

1  -  L32 

•  C  (  2 , 5  ,  X ) )  •  L3I 

80S 

C( 4 , 5 , X )  -  (C(4 , 5, X)  -  L41 

*  C(1 , 5 , X )  -  142  *  C( 2 , 5 , X } 

806 

1 

-  L43  •  C(  3 , 5 ,  X  )  )  *  L4I 

807 

C( 3 , 1  ,  X  )  -  C<  3 , 1 , X  )  -  U34 

C( 4 , 1 , X  } 

80S 

C( 2 , 1 / X )  -  C ( 2 / 1 , X )  -  C24 

C(  4 , 1 , X  ) 

809 

1  -  U23 

C( 3 , 1 , X  ) 

810 

C(l,l,f)  »  C(1,1,X)  -  U 1 4 

C  (  4 , 1  ,  X  ) 

81  I 

1  -  U13 

C ( 3 , 1 , X  )  -  U12  *  C ( 2 , 1  ,  X  ) 

812 

C( 3 , 2 , X )  «  C( 3 , 2 , X )  -  U34 

C( 4 , 2 , X  ) 

813 

C ( 2 , 2  ,  X  )  -  C( 2 , 2 , X )  -  U24 

C( 4 , 2 , X  ) 

814 

1  -  U23 

C  (  3 , 2  ,  X  ) 

815 

C( 1 , 2  ,  X)  -  C  ( 1 , 2  ,  X  )  -  U14 

C  ( 4 , 2  ,  X  ) 

816 

1  -  U13 

C( 3 , 2  , X  )  -  U12  •  C ( 2 , 2 , X } 

8)7 

C ( 3 , 3  ,  X  )  -  C( 3 , 3 , K )  -  U34 

C  (  4 , 3  ,  X  ) 

818 

C( 2 , 3 , X)  «  C( 2 , 3 , X )  -  U24 

C  (  4 , 3  ,  X  ) 

819 

1  -  V23 

C<  3 , 3 , X) 

820 

C( 1 , 3 , X )  *  C ( 1 , 3 , X )  -  U14 

C<  4 , 3 , X  } 

82  1 

1  -  U1 3 

C ( 3 , 3 , X }  -  U12  *  C ( 2 , 3 , X ) 

822 

C ( 3 / 4 , X  )  •  C( 3 , 4 , X )  -  U34 

C  (  4  .  4  ,  X  ) 

823 

C ( 2 , 4  ,  X  )  -  C( 2 , 4 , X )  -  U24 

C(4,4,D 

824 

1  -  U23 

C  (  3  ,  4  ,  X  ) 

825 

C(1 ,4  ,X)  -  C(1 ,4,X)  -  U14 

C  (  4  ,  4  ,  X  ) 

826 

1  -  Cl  3 

C( 3 , 4 , X }  -  Cl 2  *  C( 2 , 4 , X } 

827 

C  (  3 , 5  ,  X  )  -  C(  3 , 5,  X  )  -  U 34 

C  (  4  ,  5  ,  F  ) 

828 

C  7  2 . 5 , X  )  •  C ( 2 , 5 , X )  -  U24 

€  ( 4 , 5 ,  X  ) 

829 

1  -  U2  3 

C  {  3 , 5 ,  X  ) 

SOURCE  PROCRAM 

DATE  1/20/90 

PACE  0 

nse2d.f 

T1ME  1 1:4&04  am 

11 

SOURCE  TEXT 


SUBROUTINE  MATRX2 ( I  MAX , KMAX ) 


SUBROUTINE  MATKX2 


PARAMETER  <  IX-180  ,  Juc-60 ) 

C0W4ON/TR I  D1/DD<  4 , 4  ,  IX ,  EX ) 

COMMQN/TR  I D2/MM ( 4  , 4  ,  IX ,  KX } 

COHMGN/TR I D3/EE ( 4 , 4 , IX , KX ) 

COHMON/TRID4/  GG(4,IX,KX) 

COMMON/SC RAT/ A ( 4 , 4  ,  IX ) ,HB (4 , 4 , IX , XX ) ,C( 4 , 5, IX) 

REAL  MM 

REAL  U1,U1,L31,L41,L32,L32,L42,L33,L43,L44 
2  ,  L1I  ,  L2I  ,  L3I  f  L4I 

THIS  SUBROUTINE  PERFORMS  THE  BLOCK  Til  01  AGONAL  MATRIX  I  VERSION  rOR 
AH  ENTIRE  J -CONSTANT  PLANE  DURING  THE  1ETA-  SWEEP 

KM1  -  KKAX  -  1 

I Ml  -  I MAX  -  1 

DO  1  II  -1,4 

DO  1  I  -  2  ,  IM1 

AI  -  1.  /  MM(1, 1,1,1) 

GG (11,1,1)  -  GG(I1,I,1)  •  AI 

HH(I1, 1,2.1)  -  EE(I1, 1,1,1)  *  AI 
HB(  XI ,2,1,1)  -  EE (11,2, 1,1)  •  AI 
IB(X1,3 , 1,1)  -  ££(11,3,1,1)  *  AI 


HH (11,4,1 ,  1 )  -  EE ( II , 4 , I , 1 ) 
1  CONTINUE 


DO  1000  t  -  2  ,  KMAX  -  2 

DO  5  II-  1  ,  4 

DO  2  I  -  2  ,  IM1 

0(11,1,1)  -  GG( 11 , I ,K)  -  DD( II ,1,1 , K) 


-  DD( II ,1 , I , K)  •  GG(1,I,K-1) 

-  DD ( I 1 , 2 , I , K )  •  GG ( 2 , 1 , K- 1 ) 

-  DD(11 ,3,1 ,K)  •  GG ( 3 , I , K- 1 ) 

-  DD(I1,4,2,K)  •  GG( 4 , I , K-l ) 


2  CONTINUE 

DO  5  12  -  1  ,  4 

DO  5  I  -  2  ,  I  Ml 

A(I1,I2,X )-  MM( 11,12,1 , K )  -  DD ( I 1 , 1 , I , K ) 

1  -  DD ( 1 1 , 2 , 1 , X ) 

2  -  DD(I1,3,I ,K) 

3  -  DD  ( 1 1 , 4 , 1  ,  K  } 
C(ll , 12+1 , I )«  EE (II, 12, I, K) 

5  CONTINUE 


BB(1 ,12,1 , K-l) 
HH ( 2 , 1 2 , 1 , K-l ) 
BH ( 3 , 1 2 , 1 ,  K-l ) 
BH ( 4 , 1 2 , 1 , K-l ) 


DO  3  I  -  2  ,  IM1 
Lll  -  A( 1 , 1 , 1 ) 

L1I  -  1.  <■  Lll 
U12  -  A(1,2,I)  •  Lll 
U13  -  A(1 ,3,1 )  •  Lll 
U14  -  A( 1 , 4  , 1  )  *  Lll 

L21  -  A( 2 , 1 , 1 ) 

Lll  -  A( 3 , 1 , 1 ) 

L41  «  A( 4 , 1 , 1 ) 

L22  -  A( 2 , 2 , I )  -  L21  «  U12 
L21  -  1.  /  L22 

C23  -  ( A(  2  ,  3  ,  X  >  -  L21  *  U13)  *  L2I 

U24  «  ( A ( 2 , 4 , I )  -  L21  *  U14)  -  L2I 

L32  -  A( 3 , 2 , I )  -  Lll  «  U12 

L42  -  A(4 .2,1 )  -  L 41  •  C12 

L33  -  A( 3 , 3 ,1  )  -  L31  •  V13  -  L32  • 

L3I  -  1.  /  L33 

U34  «  ( A(  3 , 4 , 1 )  -  L31  *  U14  -  L32 
L43  -  A(  4 , 3 , 1  )  -  L41  •  U13  -  L42  • 
L44  -  A(4,4,I)  -  L41  *  U14  -  L42  * 
L41  =  1.  /  L44 
C(1,1,I)  -  C(1,1,I)  •  Lll 
C(  2 , 1 , 1 )  -  (C ( 2 , 1 , 1 )  -  L21  •  C(l,l 
C ( 3 , 1 , 1 )  -  (C( 3 , 1 , 1 )  -  L31  •  C(l,l 
1  -  L32  -  C ( 2 , 1 

C( 4 , 1 , 1 )  -  (C ( 4 , 1 , 1 )  -  L41  *  C(l,l 


L43  •  U34 


C  ( 1 , 2 , 1 ) 
C  ( 2 , 2 , 1 ) 
C( 3 , 2 ,1 ) 

1 

C(  4 , 2 , 1 ) 

1 

C(l,3,l) 
C(2 , 3 , I ) 
C( 3 , 3, I ) 

1 

C(  4 , 3 , 1 ) 

1 

C  ( 1 , 4 , 1  ) 
C  (  2 , 4 , 1 ) 
C( 3 , 4 , 1 ) 

1 

C  ( 4 , 4 , 1 ) 

1 

C(1 , 5, I ) 
C( 2 , 5, I ) 
C( 3 , 5, I ) 

1 

C(4 , 5,1 ) 

1 

CO, l, I) 
C( 2 , 1 , 1 ) 

1 

C(1,1,I) 

1 

C( 3 , 2 , 1 ) 
C  (  2 , 2 , 1 ) 

1 

C(  1 , 2 , 1 ) 

1 

CO,  3, 1) 
C  (  2 , 3 , 1 ) 

1 

C(  1 , 3 , X ) 

1 

C(  3 , 4 , 1 ) 
C(2, 4 , I ) 

1 

C  ( 1 , 4 , 1  ) 

1 

C  (  3  ,  5 , 1  ) 

C(2, 5,1) 


(C ( 2 , 1 , I )  -  L21  •  C(1,1.I))  •  L2I 
(C( 3 , 1 , 1 )  -  L31  •  C(1,1,I) 

-  L32  -  C (  2 , 1 , 1  )  )  •  L3I 

(C ( 4 , 1 , I )  -  L41  *  C(1,1,X)  -  L42  •  C(2,1,I) 

-  L43  •  C ( 3 , 1 , I ) )  •  L4I 

C(  1 , 2 , 1 )  *  Lll 

(C ( 2 , 2 , I )  -  L21  *  C(1 , 2 ,1 ) )  *  L2I 
(C  (  3 , 2 , 1 )  -  L31  •  C  ( 1 , 2 , 1  ) 

-  L32  *  C(2,2,I))  •  L3 I 

(C ( 4 , 2 , I )  -  L41  *  C( 1 , 2 , I )  -  L42  •  C(2,2,X) 

-  L43  *  C(3 ,2,1 ) )  •  L4I 

C( 1 , 3 , 1  )  *  Lll 

(C ( 2 , 3 , 1 )  -  L21  *  C ( 1 , 3 , 1 ) )  *  L2I 
(CO, 3, 1)  -  L31  *  C(  1 , 3 , 1 } 

-  L32  *  C(  2 , 3 ,1  )  )  •  Ul 

(C(4 , 3 , 1 )  -  L41  *  C ( 1 , 3 , 1 )  -  L42  *  C(2,3,I) 

-  L43  •  C (  3 , 3 , 1 ) )  •  L4I 


C(  1 , 4  , 1  )  •  Lll 
(C( 2 , 4 , I )  -  L21 
(C( 3 ,4 , I )  -  L31 

-  L32 

(C(4 ,4,1)  -  L41 

C ( 1 , 5 , 1 )  *  HI 
(C(2 , 5 , I )  -  L21 
(CO, 5,1)  -  L31 

-  L32 
(C ( 4 , 5 , I )  -  L41 

CO, 1,1)  *  034 

C(2,l ,1 )  -  U24 

-  U2  3 
C(1,1,I)  -  U14 

-  Ul  3 
C(3, 2, 1)  -  U3 4 
C( 2 , 2 , 1  )  -  U24 

-  U2  3 
C  ( 1 , 2 , 1  )  -  Ul  4 

-  Ul  3 
C (3 , 3 , I )  -  U34 
C(2  ,  3,1  )  -  U24 

-  U2  3 

C(1  ,3  ,J )  -  Ul 4 

-  Ul  3 
C ( 3 , 4 , 1 )  -  U34 
C(  2 , 4  , 1  )  -  U24 

-  U23 
C(  1  ,  4 , 1  )  -  C14 

-  Ul  3 
C( 3 , $ , I )  -  C34 
C{ 2 , 5 , 1 )  -  C24 


•  C  ( 1 , 4 , 1 ) )  *  UI 

•  C  ( 1 , 4 , 1 ) 

‘  C  (  2 , 4 , 1 )  )  •  L3I 

•  C ( 1 , 4 , 1  )  -  L42  •  C( 2 , 4 , 1 ) 

-  L43  •  C ( 3 , 4 , I ) )  •  L4I 

•  C(1 , 5 , I ) )  •  L2I 

•  C  ( 1 , 5 , 1 ) 

•  C(2, 5,1) )  *  L3I 

»  C ( 1 , 5 , I )  -  L42  •  C( 2 , S , I ) 

-  L43  *  C ( 3 , 5 , 1 ) )  •  L4I 
C( 4,1 ,1 ) 

C( 4 ,1 , 1 ) 

C(3,1.I) 

C(4,1,I) 

C( 3 ,1 , 1 )  -  U12  *  C( 2 , 1 , 1 ) 

C( 4 , 2 , 1 ) 

C( 4 , 2 , 1  ) 

C( 3 ,2 , 1 ) 

C( 4 ,2 , 1 ) 

C( 3,2,1  )  -  U12  *  C( 2 , 2 . 1  ) 

C  (  4 , 3 , 1 ) 

C( 4 , 3 , 1  ) 

C( 3 , 3 , 1 ) 

C(4 ,3 ,1 ) 

C( 3,3 ,1  )  -  Ul 2  •  C{ 2 , 3 , 1 ) 

C  (  4  ,4 , 1 ) 

C(4  ,4 , 1 ) 

C( 3, 4 , 1 ) 

C  {  4 , 4  , 1  ) 

C(  3 ,4 , 1 )  -  Ul 2  •  C(  2 , 4 , 1 ) 

C( 4 , 5 , 1  ) 

C(4,5.I) 


-  Ul 2  *  C( 2 , 2 . 1  ) 


C(2 , 3  » I ) 


C<  2 , 4 , 1 ) 


SOURCE  PROCRAM 

date  1/20/90 

PACE* 

nse2d.f 

TIME  11:4&04  am 

HBCJI 

SOURCE  TEXT 


SUBROUTINE  METRIC 


SUBROUTINE  METRIC 


PARAMETER  (IX-18G,k.x-60) 

COMMON/FI X/OMEGA , HDOT 

CONMON/DGR I D/DT , IMAX , XMAX , I TEL , I TEU 

COKKON/GRID1/1(IX,IX) ,Z(IX,XX) 

COMMON/GRI D/YACOB( IX , XX ) 

COMMON/MTRIX/XIX(IX , XX) ,XI2(IX,KX) , ZETAX (IX , XX) , ZETAZ( IX , XX ) , 

1XIT(1X,KX) , ZETAT( IX , IX ) 

C 

c***  THE  SUBROUTINE  METRIC  COMPUTES  THE  METRICS  IN  ALL  THE  T*0  DIRECTIONS  AND 
C  THE  QMS TE ADI  OOtmCUHTS  ETAT  ETC. 

C 

DO  1000  X  -  1  ,  IMAX 
DO  1000  1*1#  IMAX 
XTAU  -  OMEGA  *  2(1, X) 

YTAD  -  OMEGA  *  (-X<I,X))“  HDOT 
C»**  WESWT  SET  UP  IS  FOR  FLOW  PAST  AN  AIRFOIL. 


IF(I .EQ1 .OR. I  .EQ.IMAX)  GO  TO  10 
XXX  »  .5  •  (X(I41,K)-X(I-1,X)) 

ZXI  -  .5  •  (2(1+1, X)-Z(1-1,X)) 

GO  TO  IS 

10  IF(I .EQ.IMAX)  GO  TO  16 

XXI  -  1.0  •  (X(2,X)  -  X(1,X)) 

XXI  -  1.0  *  (2(2,1)  -  2(1, X)) 

GO  TO  15 

16  XXI  «  1-0  *  (X( IMAX , X )  -  X( IMAX-1 ,X> ) 

XXI  -  1.0  *  ( Z( IMAX , X )  -  Z( IMAX-1 , X ) ) 

15  CONTINUE 

IFfK.EQ.l .OR. X. EQ.IMAX)  GO  TO  17 
XZET  •  .5  *(X(X,X+1)“X(I,X-1)> 

ZZET  -  .5  • ( 2(1 ,X+1 )-2( I ,X-1 ) ) 

GO  TO  20 

17  I F  ( I  ■  EQ .  IMAX  )  GO  TO  18 

XZET  •  2.  *  X(I ,2)-l • 5  •  X(I,1)  -  .5  •  X(I,3) 

ZZET  -  2.  *  2(1,2)  -  1.5  *  1(1,1)  -  -5  •  Z(I,3) 

GO  TO  20 

18  XZET  «  1.5  •  X(I ,XMAX)-2.«  X( I , XMAX-1 )+ . 5*X ( I 
ZZET  •  1.5  •  Z(I ,KMAX)-2.*  2.( I  , XMAX-1 )+ .  5* Z< I 

20  CONTINUE 

YACOBI  •  XXI  *  ZZET  -  XZET  *  ZXI 
YACOB{ I .%)  *  1  /  YACOBI 

XIX(I,X)  -  ZZET  *  YACOB(I , K ) 

XIZ(I,X)  «  -XZET  •  Y ACOB ( I , X ) 

XTAU  -  OMEGA  •  Z(I ,  K) 

YTAU  -  -  OMEGA  •  1(1, X)-  HDOT 
XIT{ I , X )  *  -  XI X ( 1 , K )  *  XTAU  -  XIZ(I.X)  *  YTA 
ZETAX ( I , X )  «  “ZXI  •  YACOB ( I , X ) 

ZETAZ( I , X )  -  XXI  *  YACOB ( I , X ) 

ZETAT(I.K)  -  “  ZETAX (I, I)  •  XTAU  -  ZETAZ(I,X) 
1000  CONTINUE 
RETURN 
END 


X(I ,KMAX)-2.« 
Z(I ,XMAX)-2.* 


A( I , KMAX-1 )+ . 5*X ( I , XMAX- 2 ) 
1 ( I , XMAX-1 )+ . 5* Z( I , XMAX- 2  ) 


SOURCE  PROCRAM 

date  1/20/90 

nse2d.f 

T|ME  11.-4&04  am 

UNE# 


SOURCE  TEXT 


e»« 

C* 

c» 

c« 

c> 


SUBROUTINE  DISSIP 


“1082 


nwu 

TOSS 
.1086 
I-T087 
1068 

1089 

1090 

rosr 

T092 

1093 

1094 
T09S 
1096 
109? 

1098 

1099 

11 00 

no  i 

no2 

1103 

1104 

nos 
no6 
no? 

I  108 
1109 

mo 
nn 
111  2 
ni3 

1114 

ins 
If  16 
in? 
ins 

1119 

1120 
1121 
1122 

1123 

1124 

1125 

1126 

1127 

1128 
1129 

1 130 

1131 

1132 

1133 

1134 

1 1 35 

11 36 

1137 

1138 

1139 

1140 

1141 

1142 

1143 

1144 

1145 
1146 

1147 

1148 

11 49 
(ISO 

mi 

1152 

II  S3 

1154 

11 55 

1156 

1157 

1158 

1159 

1160  i  c 

1161  c 

1162 

1163  I  c 

1164  !  c 

11 65  c 

1166 

1167 

1168 
I  169 

1170 

1171 
1172! 

1173  i 

1174 

1 1 75  ; 

1176  ;  c 
I  177  c 
1178  c 


SUBROUTINE  DISSIP 


PARAMETER  (IX-180, Ax-60) 

COtMON/FLOW/QKIX.KX)  ,Q2(IX,IX)  ,Q3(IX,«)  ,0<(IX.XX) 
COf#10N/NTRIX/XIX(  IX  ,  XX  )  fXXZ(IXf  XX)  ,  XETAX(IX  ,  XX)  ,  ZETAZf  IX.  XX) . 
lXIT(IX.XX) , ZETATI IX. XX ) 

Comai/PERTR/D0I(lX,XX),D02[IX,XX).DQI(H,XX)  .DQ4(IX.XX) 
eomOW/DGRID/DT,IMAX,XMAX,ITEL.ITEU 

CaMRW/PAR/GAMMA ,  REY REF ,  ALFA ,  ALFA1 ,  REDFRE ,  AMINF ,  ALFAI 
COWCN/GR I D/Y  ACOB  ( I X  ,  XX  ) 

OOMMON/DAMP/NN ,  UNI ,  NK2X ,  NN2Y ,  NN4  X ,  m»«  V 
COMMON/RADIUS/  SPECX( IX , XX ) , SPECYIIX, XX ) 

COMON/SPEED/  U(IX,IX)  ,V(IX,XX)  .AA(IX.XX) 

DIMENSION  P(IX| .EPS (IX) , DISI ( IX, 4 ) , DIS2 ( IX , 4 ) , DISK IX. 4 1 
DIMENSION  QQ(4) 


TBIS  SUBROUTINE  ADDS  THE  FOURTH  ORDER  DISSIPATION  TERMS  TO  THE 
:  EIGHT  BAND  SIDE 


I  Ml  -  IMAX  -  1 

I M2  •  IMAX  -  2 

EMI  •  XMAX  -  1 

EM2  -  EMAX  -  2 

NDT  -  NWDT 

GM1  -  GAMMA  -  1. 

DO  20  It  *  2  ,  KMl 

COMPOTE  SWTICHING  FUNCITON  BASED  ON  SECOND  DERIVAITVE  0 F  PRESSURE 


DO  1  I  • 

P(I) 

DO  3  I  ■ 

PS1 

PS2 

EPS(I)  ■ 
CONTINUE 
EPS ( 1 ) 

EPS ( IMAX } 


1  , 

GM1 

2  ,  I  Ml 
P(I+1)  -  2 . *P ( 1  ) 
P(I  +  1  )  ♦  2. *P(I) 
ABS( PS1/PS2 ) 

-  EPS ( 2 ) 
EPS(IMl) 


IMAX 

(Q4(I, R)  -  0 . 5*Q1 { I »  X )  •  ( U  ( I ,  X )  * 


P(l-l) 


'2  +  V(I,X)**2)) 


SMOOTH  OUT  PRESSURE  COEFFICIENT 
DO  4  I  «  2  ,  IM1 

P(I)  -  AMAX1 (  EPS ( I  +  l ) , EPS ( 1 )  , EPS ( 1-1 )  ) 
4  CONTINUE 

P(l)  -  P(2) 

P(IMAX)  -  P(IM1) 


DO  10  I  -  2  ,  IM1 
C2  -  NW2X  •  P(I) 

C4  -  VW4X  -  AMIN1(»W4X,C2) 
C22  ■  C2*WDT  •  (  SPECX ( I ,  X ) 
C44  -  -C4*WDT*  SPECX ( I, X) 


SPECX ( I  +  l  ,  X )  ) 


C 

EPK 

•  04(1 

-l.E)  4  P(I-l) 

C 

HP 

*  04(1 

,X)  +  P(I) 

c 

HPP 

-  04(1 

7-1,  X)  4  P(l4l) 

DISI ( 1 , 1 ) 

-  C22' 

01(141, X)  -  Q1(I, X)  ) 

DISI ( I , 2 ) 

-  C22* 

02(1+1, X)  -  02(1, X)  ) 

DISI (I ,3) 

-  C22  * 

03(1+1, X)  -  03(1, X)  ) 

DISI ( 1 , 4 ) 

*  C22* 

04(1+1, X)  -  Q4(I,X)  ) 

c 

DISI ( 1 ,4  ) 

-  C22* 

(  HPP  -  BP  ) 

DIS2( 1 , 1 ) 

-  C44  - 

(  01(2+1, X)  -  2. *01(1, X) 

-  01(1-1, X)  ) 

D1S2(I ,2) 

»  C44  * 

(  02(1+1, X)  -  2. *Q2(I, X) 

+  02(1-1, X)  ) 

DIS2 ( I  ,3) 

-  C44  * 

(  03(1+1, X)  -  2. *Q3(I ,X) 

+  Q3 ( 1-1 , X )  ) 

DIS2 ( I  ,  4 ) 

-  C4  4  * 

(  04(1+1, X)  -  2. *Q4(I, X) 

+  04(1-1, X)  ) 

c 

DIS2(Z ,4) 

-  C44  * 

(  HPP  -  2 . *HP  +  HPM  ) 

B.C.  TREATMENT 
QO(l)  -  0K2,  X) 
Q0(2)  -  Q2<2, X) 
QQ(3)  -  03(2, K) 
00(4)  -  04(2, It) 
00(4)  -  04(2, X) 
DO  15  N4 
C2 
C22 

DISI ( 1 ,N4  ) 
DZ52(2 ,N4  ) 

IS  DIS2 ( IMAX . N4 ) 


01(1, X) 

02(1, X) 

03(1, X) 

04  ( 1  ,  X ) 

P<2)  -  04(1, X)  -  P(l) 

-1,4 

-  WW2X*P<1) 

-  C2  •  WDT  •  (  SPECX ( 1 , X ) 

-  C22  *  QQ(N4) 

-  0. 

-  0. 


+  SPECX < 2, X)  ) 


DO  16  1  -  1  ,  IM1 

DISI (1,1)  -  DISI ( I , 1 )  ^  DIS2 ( 1+ 1 , 1 )  -  DIS2 ( X , 1 ) 

DIS3 (1,2)  -  DISI ( I , 2 )  4  DIS2(I41,2)  -  DIS2(X,2) 

DIS3 (1,3)  -  DISI ( 1,3)  4  D1S2(I+1,3)  -  DIS2(I,3) 

DIS3(I,4)  -  DISI ( I ,4 )  4  DIS2 (1+1,4)  -  DIS2(I,4) 

16  CONTINUE 

1*114,  IN  DISSIPATION  TERMS 

DO  18  1  -  2  ,  IM1 

DQ1 ( I , X)  -  DIS3 ( I , 1 )  -  DI S3 ( I- 1 , 1 ) 

DO 2 ( 1 , X )  •  DIS3 ( I , 2 )  -  DIS3 ( I - 1 , 2 ) 

DO 3 ( I , X )  -  D1S3 ( I , 3 )  -  DIS3 ( I'l ,3 ) 

DQ4 ( I , X )  -  DIS3 ( I , 4 )  -  DIS3(1-1,4) 

18  CONTINUE 
20  CONTINUE 

*  DIRECTION 

DO  100  1*2,  IM1 

COMPUTE  SWCHING  FUNCTION  BASED  ON  SECOND  DERIVATIVE  OF  PRESSURE 
DO  31  X  *  2  ,  XMAX 

31  P(Jt)  -  CM1  *  (04(1, X)  -  0. 5*01(1, X)  *  (U(I,X)**2  +  V(2,X)**2)) 
P(l)  -  (  4 . *P( 2 )  -  P( 3 )  )  /  3. 

DO  33  K  •  2  ,  KM1 

PSl  -  P(X+1)  -  2 . *P ( X  }  +  P(X-l) 

PS2  «  P ( K+ 1 )  ♦  2 . *P ( X )  +  P(X-l) 

EPS(X)  -  ABS ( PS1/P52 ) 

33  CONTINUE 

EPS(l)  -  EPS ( 2 ) 

EPS (XMAX)  -  EPS(KMl) 


SMOOTH  01 T  PRESSURE  COEFFICIENT 


SOURCE  PROCRAM 

nse2d.f 


msxzii’mFi 


SOURCE  TEXT 


SmaOUTINE  SPECR  ( I  SPEC  J 


16 
r  c« 

f  s 

6  c** 

r 

T242 

1244 

1245 

Tf4« 

(247 

^9 
«  c 

T"  Q 

R 

4 


ilsr 

i* 

l1«r 
126 1' 
1262' 
1262 

1264  c 

1265  c 
1266  c 

1267  c 

1268  c 

1269  c 

1270 
1271 

(272 

1273 

1274 

1275 

1276 

1277 
(278 
1279  ! 

1280 
"1281  ' 
1282 

1283 

1284  ! 

1285 

1286  ; 
"1287  , 

(288 

1289: 

1290 ; 

(291  i 

1292 

1293 

(  294  , 

1295 

1296  ; 


SUBROUTINE  SPECS 


PMUUttTER  (II-H0,kx-«0) 

Comow/PAJ  GAMMA ,  REYREF ,  ALFA ,  ALFA1 ,  REDFRE ,  AMINF ,  ALFAI 
Q0M40N/DGR I D/DT , I MAX , KMAX , ITEL, ITEU 

C0H40W/HTR1X/  XIX(IX,KX),XIZ<IX,KX),ZETAX(IX,KX),ZETAZ(IX,KX) 

1  ,XIT(1X,IX),ZETAT(3X,KX) 

CONtON/rLOK/Ql  ( IX ,  KX )  ,Q2  ( IX ,  KX )  ,Q3  (I*  .  XX  )  ,Q4  ( IX  ,  XX ) 

C0M40N/RADI US/  SPECX<  IX  ,  KX )  , SPECY  ( IX ,  XX  ) 

C OMON/SPtED/  U(IX,KXf  , V( IX , KX) , AA(IX, XX ) 

OQMCM/GRI D/YACOB{  IX ,  XX ) 

HU  SUBROUTINE  COMPUTE  THE  SPECTRAL  RAD r US  FOR  SCALING  TEE 
XXRXJC1T  AND  IMPLICIT  DISSIPATIONS 

tfCL  -  I  MAX  -  1 

KM1  -  KMAX  -  1 

CAM  -  GAMMA  *  (  GAMMA  -  1.  ) 

DO  5  K  -  1  ,  KMAX 

DO  5  I  -  1  ,  IMAX 

0(1, K)  -  Q2(l,l}  /  Ql(l.K) 

V(I ,K)  •  03(1. K)  /  Q1 < I , K ) 

AA(I  ,f  )•  CAM  *  (  Q4  ( I ,  K )  ,/Ql  ( I ,  K )  -  0. 5*  (0(1 ,1)^*2  +  V(I,K)**2)  ) 
IF(AA( I , K ) . LT. 0 . )  PRINT* , 'NEGATIVE  A*A  -  * , AA(I , K) , 1  AT  *,I,K 
5  CONTINUE 


COMPUTE  IMPLICIT  DISSIPATION  SCALING 
SPECX  -  SPECTRAL  RADIUS  FOR  XI-DIRECTION 
SPECY  -  SPECTRAL  RADIUS  FOR  ZET A -DIRECT! ON 


IP( ISPEC . EQ. 1 
DO  20  F 
DO  20  1 

SPECX ( I , K  J 
SPECY (I, K) 
20  CONTINUE 
ELSEir( ISPEC. 
DO  30  K 
DO  30  I 
UCON 

VCON  • 

SPECX ( I ,  K  ) 
SPECY ( I , K ) 
30  CONTINUE 
ELSEIF( ISPEC. 
DO  40  K 
DO  40  I 
UCON 
VCON 
XI 2 
ZETA2 

SPECX(I.F) 
SPECY ( I , K ) 
40  CONTINUE 
ENDIF 
RETURN- 
END 


)  THEN 

-  1  ,  KMAX 

-  1  ,  IMAX 

«=  1  .  /  YACOB(I.K) 

-  1.  /  YACOB(I.K) 

EQ . 2 )  THEN 

-  1  ,  KMAX 

«  1  ,  IMAX 

U(I,K)*XIX(I,K)  ♦  V(I , K  )«XIZ(I , K)  4  XIT(I,K) 

U(I  ,K)*ZETAX(I  ,K)  4  V(X,K)*ZETAZ(X,K)  4  ZETAT(I.K) 

-  ABS(UCON)  /  Y ACOB { I  ,  K ) 

-  ABS(VCON)  /  YACOB(I.K) 

EQ.3)  THEN 
«  1  ,  KMAX 

-  1  ,  IMAX 

-  V( I , K ) *XIX ( I , K )  *  V(I,K)*XIZ(I  ,K)  4  XIT(I.K) 

-  U ( I , K ) *  ZETAX ( I , K )  +  V ( I , K ) • ZETAZ ( I , K )  4  ZETAT ( I , K ) 

-  XIX ( I , K ) *  *  2  4  XIZ(I,K)**2 

-  ZETAX  (I  »K)**2  -  ZETAZ(  I ,  K  )  •  •  2 

-  (ABS(UCON)  4  SQRT(AA(I,K)*XI2>  )  /  YACOB(I.F) 

=  (ABS(VCON)  4  SQRT ( AA( I , K) * ZETA2 )  )  /  YACOB(I.K) 


SOURCE  PROGRAM 


DATE  1/20/90  PACE* 

"M*  1 1:4&04  am  17 


SOURCE  TEXT 


J2JJ  C*« 
1249  c* 
300  c* 
“1101  c* 
1102  c»> 
1303 
*304 
1305 
T306 

1107 

1108 

1109 

1110 
in  i 
1112 

Ifri 

m 

■nr  r 

TIT* 

1119 

1120 

1321 

1322 
1123  c 

T  324 
132$ 

1326 
1127  c 
1 32*  c 

1329  C 

1330 

1331 

1332 

1333 

1334 
1 3  3  S 

1336 

1337  ! 

1338 

1339 
13  40  ; 

1341 

1342 

1343 

I  344  •  C 
134$ 

1346 
I  347  ] 

I  348  I 
1349 
i3so ; 

1 3  S I  1 
1  352  ; 

1353 

1354 

1355 

1356 

1357  ; 

1358 

1359 

1360  . 

1361 

1362 

1363 

I  364  , 

1365  ; 

1366 

1367 

1368 

1 369  i 

1370  c 

1371  c 

1372  c 
(373 

1374 

1375 
1376  c 

(177  c 
137*  c 
1179  c 
13*0 
13*1 
13*2  c 
13*3 
13*4 

13*6  ! 

1 3*7  i  c 
1388  c 
13*9  c 

1190 

1191  i 

1392  j 

1393 

1194  j 

1195 

1196  , 

1397  ; 

1398 

1399  ' 

1400 

'  1401  ,  C 

1402  c 

1403  C 

1404  :  c 

1405 

1406 
1407 

1408  j  C 

1409  C 

1410  ;  c 

1411  C 

1412  ;  C 

1413  ;  C 

1414  j  C 

1415 

1416 


SUBROUTINE  tlPBC(Cl) 


SUBROUTINE  EXPBC 


PARAMETER  ( IX-180 ,  Xx-tO ) 

COKHCN/SURF/PSUR  (II) 

C0KN0N/GRID1/X( IX , XX ) , Z(IX, XX) 

COMMON /P AP/GAMMA , REYREF , ALFA , ALFA1 , REDFRE ,  AMI  N  F 
COMHON/DGRI D/DT , IMAX , KKAX , I TEL , I TEU 
CO6M0N/CR  I D/Y ACOB  ( I X ,  IX ) 

COMNON/DAMP/WW ,  HWI  ,  WW2X  ,  WW2Y  ,  WW4X  ,  WW4Y 

CONNON/HTRIX/XIX ( IX ,  XX )  ,XIZ(  IX  ,  XX )  ,  ZETAX(  IX ,  XX)  ,  ZETA2  ( IX ,  XX) , 
1XIT( IX ,  XX  )  ,  ZETAT(  IX  r  KX  ) 

COMMON/rLOK/Ql(IX,KX),Q2(IX,KX)  ,Q3(IX,KX)  ,04  (IX,  XX) 
COMMON/SPEED/  U< IX , KX ) , V( IX , XX ) , AA( IX , XX ) 

COMMON/IN1TI/UINF , VINF 
COMMON /BC LOG/C I RCOR 
COMMON/LCK.IC/RSTRT ,  PITCH ,  RAMP 
LOGICAL  RSTRT, PITCH, RAMP 
LOGICAL  CIRCOR 
DIMENSION  Cl ( 4 ) 

DIMENSION  A(2,2) ,RHS<2) 

DATA  PI/3 . 1415927/ 

DATA  CHORD /I. / 

GAMI  “  1 .  /  GAMMA 
CAMM  •  GAMMA  -  1. 

GAMKJ  -  X.  /  GAMM 

INVISCID  AND  VISCODS  B.C.  ON  SOLID  WALL 

DO  1  I-  ITEL  ,  1TEU 

X  •  3 

Cl(l)  *  XIT( 1 , X ) 

Cl ( 2 )  -  XIX(I.X) 

Cl ( 3 )  -  XXZ(lfX) 

UCON3  *  <Q2(1,X)«C1<2)+Q3(I,X)*C1(3)) 

1/01(1, X) 

X  -  2 

Cl(l)  -  XIT( I , X) 

Cl ( 2 )  -  XIX ( I , K ) 

Cl  ( 3 )  -  XI Z ( I , K ) 

UCON2  -  (Q2 (I ,  X)  *C1 ( 2 )+Q3 (I , X)*C1(3) ) 

1/01 ( I , X) 

RHS(l)  •  2.  •  UCON2  -  UCON3  -  XIT(I,1) 

FOR  VISCOLS  FLOWS  SET  UCON  TO  ZERO  ALSO 
IFfREYREF  GT.O.)  RHS(l)  -  -  XIT(I,1) 

A( 1 , 1 )  -  XIX(1 ,1) 

A( 1 , 2 )  -  XIZ(I ,1) 

A( 2 , 1 )  -  ZETAX (1,1) 

A(2 ,2)  -  ZETAZ (2,1) 

RHS ( 2 )  -  -  ZETAT (1,1) 

TEMPI  -  A( 1 , 1 ) 

TEMP 2  -  A( 1 , 2 ) 

TEMP  3  -  A  (  2 , 1 ) 

TEMP 4  -  A( 2 , 2 ) 

DEN  -  1.  /(TEMPI  *  TEMP 4  -  TEMP 2  •  TEMP 3 ) 

A(  1 , 1 )  •  A(  2 , 2  )  •  DEN 
A( 1 , 2 )  -  -  TEMP 2  •  DEN 
A( 2 , 1 )  -  -  TEMP 3  •  DEN 
A(  2 , 2  )  -  TEMPI  •  DEN 


01(1,1)  »  2.  •  01(1,2)  -  01(1,3) 

Q2 ( I  ,  1  )  -  Q1(I,1)*(A(1,1)*RBS(1)4A(1,2)*RHS(2)) 

03(1,1)  -  Q1(I,1)  *  ( A(  2 , 1 )  *RHS  ( 1 )  •*  A(  2 , 2  )  *RHS  (  2  )  ) 

1  CONTINUE 

DO  10  I* I  TEL  , ITEU 

P2-GAMM* ( 04 ( 1 , 2 ) -  0 . 5*01 ( 1 , 2 ) • (U ( 1 , 2 ) * »2+V( 1 , 2 ) * • 2 ) ) 

P3-GAMN*  (  04  ( I  ,  3  )  -  0 . 5*01  ( I  ,  3  )  •  (  U  ( 1 , 3  )  *  *  2-9  V  ( 1 , 3  )  •  *  2  ) ) 

PI* (4 . *P2“P3 )/3 . 

PSUR( I )" ( GAMMA* PI -1 . )/( .7*AMINF**2) 

10  Q4(I ,1 )-Pl/GAMMfO. 5*Q1(I ,1)* (U(I , 1)**2+V(I , 1 ) » • 2 ) 

FAR  FIELD  BOUNDARY  CONDITION  ONLY  FOR  STEADY  FLOW 

CIRC  -  0 

ir< PITCH. OR. RAMP)  GO  TO  999 
IF(AMINF.GT.l. )  GO  TO  65 

CIRCULATION  CORRECTION  AT  THE  FAR  FIELD  IS  BASED  ON  POTENTIAL 
VORTEX 

BETA  -  SORT(  1.  -  AMINF**2  ) 

IP(CIFCOR)  CIRC  -  0.25  *  CHORD  •  CL  *  BETA  «  AKINF  /  PI 

C05AL  -  COS (ALFA) 

SINAL  -  SIN (ALFA) 

AINF  -  1. 

BINF  -  GAMMI  +  0.5  *  AMINF*«2 

CIRCQ  -  CIRCULATION  CORRECTION  QUANTITY 

X  -  XMAX 

DO  60  I  -  2  ,  IMAX- 1 

XLOC  -  X(I,X)  -  XREF 
ZLOC  -  Z(I,X) 

RADIUS  *  SORT (  XLOC*  *2  ZLOC*  *2  ) 

ANGLE  -  AT AN2 ( ZLOC, XLOC) 

CIRCQ  •  CIRC  /  (  RADIUS  *  (  1.  -  ( AMINF*SIN ( ANGLE- ALFA ))♦ *2  )  ) 
UF  •  UINF  +  CIRCQ  •  SIN ( ANGLE ) 

VF  *  VINF  -  CIRCQ  *  COS (ANGLE) 

AF5Q  *  CAMM  *  (  HINF  -  0.5*(  UF**2  +  VF**2  )  ) 

Ar  •  SQRT(AFSO) 

NOHREFLRCTING  B.C.  BASE  ON  1~D  RIEHANN  INVARIANTS 
ZTTXN,  IETZN  *  NORMALIZED  NORMAL  VECTOR 

ANOR  -  1.  /  SORT (  ZETAX (2 , X) • *2  *  ZETAZ ( I , X) * *2  ) 

ZETXN  -  ZETAX ( I, X)  *  ANOR 
ZETZN  -  ZETAZ ( I, X)  *  ANOR 

CBECX  FOR  INFLOW  OR  OUTTLOW 

FOR  INFIOK-  *1,  VEIT.  ENTROY  ARE  SPECIFIED  AS  FREE  STREAM  VALUES 
R2  IS  EXTRAPOLATED 

FOR  OUTFLOW;  R1  15  SPECIFIED  AS  FREE  STREAM  VALUE 
R2,  VELT,  ENTROPY  ARE  EXTRAPOLATED 

RHOEXT  -  01(1,1-1 ) 

RBO!  *  3  ,  01 ( J,x-i ) 


SOURCE  PROGRAM 

nse2d.f 


2SIE 


UNE  * 


SOURCE  TEXT 


DAT!  1/20/90 

PACE* 

time  1 1:46*4  am 

18 

1417  ] 
141ft  ! 
1419  1 
WO  ; 
U2I  c 

1422  c 

1423  !  C 

1424  [  c 

1425  ! 

1 42ft  i 

1427 

1428 

1429 

1430 

1431 
U32 
141$ 

1434 

1435 

1436 

1437 
143ft 

1439 

1440 

1441 
T442 

1443 

1444 

1445 

1 446 

1447 
144  ft 

1449 

1450 


UEXT  -  '■:(  1,1-1)  •  RHOI 

VEXT  *-  J3(I,R-1)  *  RHOI 

EEXT  -  Q4(I,R-i) 

PEXT  -  GAMM  *  (  EEXT  -  0.5*RHOEXT*(  UEXT* *2  +  VEXT**2  )  ) 

SET  RIEMANN  INVARIANTS  tl,  AND  R2 

VEIN  -  NORMAL  VELOCITY,  VEIT  -  TANGENTIAL  VELOCITY 

R1  -  ZETXN  •  Ur  4  ZETZN  •  VF  -  2 .  *  AF  *  GAMMI 

R2  »  ZETXN  *  UEXT  4  ZETZN  •  VEXT  42.*  SQRT (  GAMMA  •  PEXT  / 

1  RHOEXT  )  *  GAMMI 

VELN  -  (  R1  4  R2  )  *  0.5 

SPSQ  -  (  R2  -  R1  )  *  GAMM  «  0.25 

A2  -  SPSQ*  *2 

SET  OTHCR  FIXED  OR  EXTRAPOLATED  VARIABLES 

IF(VELN.LE.O. )  THEN 

VELT  -  ZETZN  *  UF  -  ZETXN  •  VF 
ENTRO  -  GAMMA 
ELSE 

VELT  -  ZETZN  *  UEXT  -  2ETXN  *  VEXT 
ENTRO  -  RHOEXT*  * GAMMA  /  PEXT 
ENDir 

NON  COMPUTE  FLOW  VARIABLES 


'  ! 


U(J,I)  - 
V(X,R)  • 
01 (I, R)  ■ 
PRESS 
Q2< I i K)  * 


2ETXN 

ZETZN 


1451 

Q4 ( 1 , R ) 

1  452 

60 

CONTINUE 

1453 

GO  TO  67 

1 454 

65 

CONTINUE 

U55 

C 

1456 

C 

B.C.  FOR 

1457 

C 

1 458 

R  -  RMAX 

1459 

DO  66 

1460 

RHOI  - 

VELN 
VFL.S  - 
(  A2  «  ENTRO 
A2  •  01(1 /R) 
Ol(IpK)  *  0(1, R) 
Q1(I,R)  *  V( I , R ) 
PRESS’GAMMI  >  0.5 


ZETZN 

ZETXN 

GAMI 

GAMI 


VELT 
VELT 
•*  GAMMI 


QI ( I , R  }  •  (  0(1, R)* 


♦  V(I ,R)**2  ) 


2  ,  IMAX-1 

01(1 -R) 


1461 

U(I ,R) 

Q2 ( 1 <  K )  *  RHOI 

146? 

V(I ,R) 

Q3(I,R)  •  RHOI 

1463 

ANOR 

1 .  /  SQRT (  ZETAX ( 1 , R ) * 

l  464 

ZETXN 

ZETAX ( I , R )  *  ANOR 

1465 

ZETZN 

ZETAZ(I.K)  •  ANOR 

1466 

VELN 

ZETXN  •  U( I , R)  4  ZETZN 

(467  - 


IF(VELN.CE.  0.  )  THEN 


•2  4  ZETAZ( I , R ) * • 2  ) 


V(I,R) 


1468 

Ql(I,R)  -  Ql(I.R-l) 

1469 

Q2(I,R)  -  Q2(I,R-1) 

1470 

Q3 ( I , K )  -  Q3 ( I , R- 1 ) 

1471  . 

Q4fI,R)  *  Q4  ( 1  ,  R-*I ) 

1472  ; 

ENDIF 

1 473 

66 

CONTINUE 

1474  c 

1475 ; 

67 

CONTINUE 

1476  c 

1477  c 

OUTFLOW  B.C.  AT  THE  DOWNSTREAM  OF  C-GRID  ONLY  FOR  INVISCID 

i  1478  C 

i  1 479  , 

IF(R£YREF . LT. 0. )  THEN 

1480, 

I  -  1 

,  1481 

IP1  -  1 

1  148? 

SIGN  *  -1 . 

1483  ;  c 

;  1484  j  c 

CHECR  FOR  SUPERSONIC  FLOW 

1 485  •  C 

I486 

I F( AMINF . GT . 1 . )  GO  TO  75 

i 487  ;  C 

I486 

72 

CONTINUE 

1489  i 

DO  74  R  -  1  ,  RMAX 

1 490  ;  c 

'  1491  1  c 

CORRECT  FREE  STREAM  VALUES  WITH  CIHCULATION  CORRECTION 

1492  j  c 

1493 

XLOC  -  X(I , R )  -  XREF 

1494 

ZLOC  -  Z(I ,R) 

1 495  | 

RADIUS  -  SQRT (  XLOC**2  •»  ZLOC**2  ) 

1496  ; 

ANGLE  -  ATAN2( ZLOC, XLOC) 

1497 

1498 

1 499 

1500 

1501 

1502 

1503 

1504 

1505 
\  50ft 
1507 

1 50ft  i 
1509  ! 
15)0  ! 
1511  ■ 
I  St  2 
1513  j 
I  S 1 4 

1515  I 

1516 

1517 ; 

1516  , 

1519  1 

1520  ; 

1521 

1522 

1523 

1524 

1525 

1526 

1527 

1528 

1529 

1530 

1531 

1532 

1533 

I'  14 

’  35 
I  >36 


CIRCQ  -  CIRC  /  (  RADIUS  •  (  1.  -  ( AMINF*SIN ( ANGLE-ALFA) ) • *2  )  ) 
UF  •  UINF  4  CIRCQ  *  SIN (ANGLE) 

VF  -  VINF  -  CIRCQ  *  COS (ANGLE) 

AFSQ  -  GAMM  *  (  HINF  -  0.5*(  UF**2  4  VF**2  )  ) 

Ar  -  SQRT (AFSQ) 

XXXB,  XI  ZH  -  NORMALIZED  HORIZONTAL  VECTOR 

ANOR  -  1.  /  SOFT (  XIX(I,R)**2  4  XIZ(I,R)**2  ) 

XIXH  -  X1X(I , R )  *  ANOR 
XI ZH  -  XI Z ( I , R)  »  ANOR 

RHOEXT  •  Q1(I4IP1,R) 

RHOI 
UEXT 
VEXT 
EEXT 
PEXT 


-  1 .  /  Ql(I*IPl,K) 

-  Q2(I4IP1,R) 

-  Q3 ( 1 4XP1  ,  R ) 

-  Q4 ( I +IP1 , R ) 

-  GAMM  *  (  EEXT  -  0.5*REOEXT*(  UEXT* 


2  4  VEXT*  *  2  )  ) 


SET  RSEMANN  INVARIANTS  R1 ,  AND  R2 


R1  -  XI1H  *  OF  ♦  XIZH  •  VF  -  SIGN  *  2. 
R2  -  XI IB  •  CEXT  *  XI ZH  *  VEXT  *  SIGN 


•  AF  *  GAMMI 

2  •  SQRT (  GAMMA  •  PEXT  / 

RHOEXT  )  •  GAMMI 


VELN  -  (  FI  ♦  R2  ) 
SPSQ  -  SIGN  •  (  R2 
A2  -  SPSQ* *2 


0.5 
FI  ) 


GAMM  •  0.25 

SET  OTHER  FIXED  OR  EXTRAPOLATED  VARIABLES 


I F( SIGN* VELN . LE . 0 . )  THEN 

VELT  -  XI  ZB  *  UF  *  XIXH  •  VF 
ENTRO  -  GAMMA 
ELSF 

VELT  •  - XI ZH  *  UEXT  ♦  XIXH  •  VEXT 
ENTRO  -  RHOEXT* *GAMMA  /  PEXT 
ENDIF 


C 

C 


COMPOTE  FLO*  VARIABLES 


SOURCE  PROGRAM 

nse2d.f 


date  1/20/90  PACE  * 
TIME  ll:4&04*m 


LINE# 


1537  i  c 
“1536 
1539 
“1540  ' 

I54|  ! 

1542  s 

1543  : 
1544 ; 
1545 


SOURCE  TEXT 


XIXH  *  VELN  -  XIZH 
XI2.H  •  VELN  +  XI XB 


Q1  ( I , K )  »  (  A2  *  ENTRO  •  GAMI  )  **  GAKM1 
PRESS  -  A2  *  01(1,*)  »  GAMI 
02  ( I ,  K )  -  Ql(I.K)  *  0(1, K) 

03(1,*)  -  01{I, *)  •  V( 1 , K ) 

Q4(I.*)  -  PRESS*GAMM1  4  0.5  «  01(1,*)  * 
74  CONTINUE 


<  0(1 ,  *  )  *  *  2  ♦  V( I , * ) «  *  2  ) 


1546  1 

IF(I .EQ. IMAX)  GO  TO  79 

1547  1 

I  -  IMAX 

1548 

IP1  -  -1 

1549 

SIGN  -  1. 

1550 

GO  TO  72 

ISSI 

79  CONTINUE 

1  5S2 

GO  TO  89 

ISS3  C 

>5.54  c 

B.C.  FOR  SUPERSONIC  FLOW 

TSSS  c 

1  5S6 

75  CONTINUE 

1 557 

DO  84  I  -  1  ,  KMAX 

I5S8  C 

1559 

RIO I  -  1.  /  OKI,*) 

1  560 

0(1,*)  -  02(1,*)  *  RHOl 

1561 

V(I ,*)  -  03(1 ,*)  *  RHOI 

1562 

AMOR  -  1  /  SQRT(  XIX(1 

1563 

XIXH  -  XIX(I , K )  •  ANOR 

1564 

XI ZB  -  XIZ<I ,*)  •  ANOR 

1565  i 

VELN  -  SIGN  •  (  XIXH  • 

1566  1  c 

1  567 

I F( VELN . GE ■ 0 . )  THEN 

t  568 

01(1 , K )  «  Q1 ( I ♦ IP1 , K  ) 

1569 

Q2 ( I , * )  *  02(I*IP1 ,*) 

1570 

03(1,*)  -  Q3(J4lPl,K) 

1571  ’ 

Q4(l  ,*)  -  Q4 ( I  *  I  PI , *  ) 

1572 

ESDI  F 

1573 

84  CONTINUE 

1574 

1F(I .EQ. IMAX )  GO  TO  89 

1575 

I  -  IMAX 

1576 

IP1  -  -1 . 

1577 

SIGN  •  1 

1578 

GO  TO  75 

1579 

89  CONTINUE 

0(1 ,*)  ♦  XI ZB  •  V(I ,*)  ) 


1580  c 

1581  C 

1582  c 

1583 

1584  \  c 

1585  C 

1586  C 

1587  •  C 

1588  C 

1589  ' 

1590 

1591 

1592  ' 

I  593 

1594 

1595 

1596 

1597 

1598 

1599 

1600 
1601 
1602 

1603  ‘ 

1604  , 

1605  C 

1606 
1607 ;  c 

1608  ,  C 

1609  C 


DOWNSTREAM  B.C.  FOR  VISCOUS  FLOW 
EISEI F(RFYREF .  GT .  0.  )  THEN 
OUTFLOW  B.C.  FOR  VISCOUS  FLOW 

RBO ,  V  AND  V  ARE  SET  TO  THE  VALUES  OF  THE  NEXT  INTERIOR  POINTS 
PRESSURE  IS  EXTRAPOLATED  AND  THEN  COMPUTE  ENERGY  <Q4) 

DO  100  *  «=  1  ,  KKAX 

I  -  1 

old,*)  -  oi(i4i,*> 

02(1,*)  -  Q2(l*l,*) 

03(i,*)  -  03(i*i, rj 

PEXT  •  GAKM  •  (  Q4 ( 1  +  1 , * )  -  0  5  *  (02(1*1  ,K)**2  4  Q3 ( 1*1 , * ) *  *  2  ) 

1  /  01(1*1,*)  ) 

04(1,*)  -  PEXT.  GAMM  *  C  5  *  <  02(1, X)"2  *  0J(X,*>««2  )  /  01(1,*) 

I  «  IKAX 

Q1(I,*)  -  OKI-1.  *) 

02(1, *>  -  02(1-1,*) 

Q3 (I , * )  •  03(1-1,*) 

PEXT  -  GAMM  •  (  04(1-1,*)  -  0.5  •  (02 ( I - 1 , * ) • * 2  *  0*( 1-1 . K ) • *2 ) 
1  /  01(1-1,*)  ) 

04(1,*)  -  PEXT /GAMM.  4  0.5  •  (  02(1,*)*  *2  *  Q3(I,K)**2  )  /  Ql(l,*) 
100  CONTINUE 


AVERAGE  FLOW  VARIABLES  ACROSS  WARE  CUT  (FOR  C-GRID) 


99  9  CONTINUE 

DO  120  I  -  1  ,  I  TEL  -  1 

IU  *  I MAX  4  j  -  i 

Q1AVG  -  05  •  (  Q1(I ,2)  4  01 ( IU , 2 )  ) 

Q2AVG  -  0  5  •  (  Q2(1.2)  *  02(11,2)  ) 

03AVG  -  0.5  •  (  01(1.2)  4  Q3 ( IU , 2 )  ) 

01(1,1)  ■  01AVG 

01(IU,1)  ■  Q1AVG 
02(1,1)  -  02AVC 

02 (It'd)  •  Q2AVG 
03(1.1)  -  Q3AVG 

Q3 ( IU , 1 )  •  Q3AVG 

PSLON  •  GAMM  *  (  04(1,2)  -  0.5  *  < 


04(1,1)  -  Q4AVG 

Q4 ( I U , 1 )  -  04AVG 
120  CONTINUE 
RETURN 
END 


GAMM  »  (  04(1,2)  -  0.5  *  (  Q2(I,2)*«2  4  03(I,2)**2  j  / 
01(1.2)  ) 

GAMM  *  (  04 ( I U , 2 )  -  0.5  •  (  02 ( Tf . 2 )  • 24Q3 ( IU , 2 ) • • 2  )  / 
01(JU,2)  ) 

0  5  *  (  PS  LOW  ♦  PSUPP  ) 

PSAVG  /  GAMM  *  0.5  «  (  02(I,1)«*2  4  Q3(I,1)«*2  )  / 

Qld.l) 


SOURCE  PROCRAM 

nse2d.f 


SOURCE  TEXT 


SUBROUTINE  STRESS ( ITN ) 


•  637  j  c« 

1638  '  C- 

1639  |  C” 


SUBROUTINE  STRESS 


I 

1 64  S  ' 

1649  ! 

1650  1 

1651  c 

1652  c 

1653  c 

1654  c 

1655 
1656 ! c 
(657  ; 

1658 

1659 

1660  I  c 

I  tool  C 

1662  !  c 

1663 

1664  ' 

1665  i 

1666 

1667 

1668 
1669 
16/0 

1671 

1672 

1673  ' 


PARAMETER  < IX- n 0 , kx *60 ) 

C0M10N/PL0H/Q1  (XX,  KZ)  ,Q2(IX,FX)  ,QJ(IX,KX),Q4(  21,  «) 
COWfON/DGRID/DT  ,  IMAX  ,  KMAX ,  I  TEL ,  ITEU 
COtMON/GRLDl/X  ( IX  ,  KX  )  , 2( IX , KX ) 

COIPION /PAR /GAMMA,  REYREF,  ALFA ,  ALFA1  ,REDFRE,  AMI  NF,ALFAI 
COMMON /SPEED/  V ( IX , XX ) , V( TX , XX ) , AA{ IX , KX ) 

COM-ON/PER TR/D01 ( IX, XX ) ,DQ2< IX, XX) ,DQ3(IX,XX) ,DQ4(*X,IX) 
C06040N/KUTUR/CMU  ( I X  ,  XX  ) 

DIMENSION  RBI (IX), RH2 (IX),RB3(IX) ,RH4(IX) 

COMMON/ LOG I C/RSTKT , PITCH , RAMP 
LOGICAL  RSTRT, PITCH, RAMP 

Til 6  SUBROUTINE  ADDS  VISCOUS  TERMS  TO  THE  BIGHT  HAND  SIDE 


CMDL  -  LAMINAR  OR  MOLECULAR  VISCOSITY 

CMUL  -  AMI.  -  /  REYREF 
CALL  EDDY (CMUL) 

KMAXM1  -  XMAX  -  1 
IMAXM1  -  IMAX  -  1 

PR  -  1. 

COMPUTE  TXX , TXY  AND  VISCOUS  DISSIPATION  AT  I  -  1  /  2 

DO  30  X  -  2  ,  XMAXM1 
DO  20  I  •  2  ,  IMAX 
UXI  -  U(I ,X)  -  Ltl-l.X) 

VXI  •  V(1 ,X)  -  V<1-1,X) 

AXI  -  AA ( 1 , X )  -  AA<I-1 ,X) 

UZET-  .25*(L'(I,X4l)-V(I,X  !’ '  1-1 ,  X*1  )-U(  1-1 ,  X-l )  ) 

VZET-  .  25 * { V( I  ,X*I)-V(I  I-1,X-H)-V(I-1  ,X-1)  ) 

AZET-  .25*  (  AA(  1  ,X-*1)~AA<I  ,  X-l ) +AA(  1-1 ,  X+l  )- AA(  I-I ,  X- 1  )  ) 
XXI  -  X(I .X)  -  X(I-l.X) 

ZXI  -  Z<I ,X)  -  1(1-1, X) 

XZET-  .25  •  (X(I,K-*1)-X(I,X-1)-»X(I-1,X+I)-X(I-1,X-1)) 


1674 

ZZET- 

.25 

•  (Z(I  ,F-*1)-Z(I  ,K-1) 

1675 

YAC 

■ 

1X1 

•  ZZET  *  ZXI  • 

XZET 

1676 

YAC 

■ 

1 

YAC 

1677 

XIX 

■ 

ZZLT 

•  YAC 

1678 

ZETAX- 

-  ; 

XI  *  YAC 

1679 

XIZ 

• 

-x:et  -  YAC 

1680 

ZETAZ- 

XXI 

*  YAC 

1681 

CNM 

(  5 

•  (CMV(I.F)  4 

CMV  ( I 

1682 

ux 

UXI 

•  XIX  -  UZET  • 

ZETAX 

1683 

VX 

VXJ 

•  XIX  *  VZET  • 

ZETAX 

1684 

AX 

AXI 

•  XIX  4  AZET  ♦ 

ZETAX 

1685 

UZ 

UXI 

•  xi :  4  tzet  • 

ZETAZ 

1686 

vz 

VXI 

•  XIZ  *  ET  • 

ZETAZ 

1687 

A2 

AXI 

•  XIZ  4  AZET  • 

ZETAZ 

1688 

TXX 

-(-4 

.  •  UX  *  2  • 

VZ)  • 

1689 

TXY 

CNM 

•  (UZ  *  VXI 

1690 

TYY 

“CN  v 

/  J  4  (-4.  • 

V7  4 

1691  C 

R4 

m ' 

X,f  )«U(J-1,X)) 

•TXX*  ( 

1692  ;  c 

1 

CN  v 

/  PR/ (GAMMA  - 

1)  * 

1693  C 

1694  c 

1 695  c 

1 696  c 
1  697 

1698 

1 699 

1700  . 

1701  ; 
'702 
1703 
1704 ; 

1705 

1706 

1707 

1 708  ; 

1 709  j  c 

1710 

1 7 11  | 

1712  , 

1713 

1714  ! 

1715 

1716 

1717 

1718 

1719 

1720 

1721 

1722  j 

1723  , 

1724  j 

1725  j 

1726  | 

1727  1 

1728  ] 

1729 

1730  | 

1731  | 

1732  i 

1733 

1734 

1735  I 

1736 

1737  j 

1738  |  c 

1739  c 

1740  c 

1741  j  c 

1742 

1743 

1744 

1745 

1746 

1747 

1748 

1749 

1750 

1751 

1752 

1753 


,  XMAX 

(C(I+1  ,X)-U<I-l  ,K)+ti<I  +  l  ,X-l)-U<I-l,X-l)> 

( AA(  I  ♦  1 ,  X  ) -AA(  I  - 1 ,  X  )■*  AA  ( I  *  1  .  X- 1  >-  AA<  I- 1 ,  X- 1  ) ) 
<X<I*1,X)-X(I-1,K)*X(I+1,K-1)-X<1-1,X-1)) 
<Z(I+1  R)-Z(I-1,K)+Z(1*1,K-1)-Z(X-1,K-1)) 

)  -  ;<i,x-i> 


54  -  (<!  (I  ,X)*U<J-i  ,X)  )*TXYMV(I,X)+V<3-l,X))*Tnr)*0.5 

1  -  CNV  /  PR  /  (GAMMA  -  1.)  *  AZ 

DEBUG 

TURN  OFF  ENRGY  DISSIPATION  AND  DIFFUSION 
R4  -  0 
54  -  0 
BB 1(1)  -  0. 

RB2 ( I  )  *  (XIX  *  TXX  *  XIZ  •  TXY)  /  YAC 
FH 3(1)  «  (XIX  •  TXY  -  XII  »  TYY )  /  YAC 
20  RH4 ( I )  ■  (XIX  •  R4  -  XII  *  S4)  /  YAC 
DO  JO  I  -  2  ,  1MAXMI 

DQ1 ( I , X)  -  DQl(I.X)  *  RBI (1*1)  -  RB1(I) 

DQ2  ( I  ,  X  )  -  DO  2  ( I  ,  X  )  •»  RH2(I-*1)  -  RB2  <  2  ) 

DQ3 ( I , X )  -  DQJ(X.X)  -  RlJ(X-l)  -  RB3(I) 

DQ4  ( I  ,  X )  -  DQ4  (  2  ,  X  )  ■»  R|4(I*1>  -  FB4(J) 

JO  CONTINUE 

IN  THE  l  DIRECTION 
DO  70  I  -  2  ,  IKAXM1 
DO  60  X  -  2  ,  XMAX 

UXI  -  .25  •  (U(I*1 ,X)-U<I-1 ,K)+(H1+1 ,*-l)-C(I-l,lt-l)> 

VXI  -  .25  *  (V(I*U)-V(I-l,r)-V(W(M)-V(l'l,M)) 

AXI  -  .25  •  (AA(I-»1,X)-AA(I-1,X)  +  AA(I-»1.X-1)-AA(I-1,X-1) 
XXI  -  .25  •  <X(Z-»2,K)-X(I'1,X)+X(7-»1,X-:L)-X<1-1,E-1}) 

ZXI  -  .25  •  (Z(I  +  1  R)-Z(I“1  ,X)  +  Z(I-*1,X-I)-Z(I-1,X-1)) 
DIET  -  U(I ,X)  -  V< 1,1-1) 

VZET  -  V(I ,X)  -  V(X,X-1) 

AZET  *  AA ( 2  »  X )  -  AA(7,K-2) 

XZET  -  X( I , X )  -  X(I ,X-1) 

ZZET  -  If  I ,X)  -  1(1 , X-l ) 

YAC  ■  XXI  •  ZZET  -  171  •  XZET 

YAC  -  l .  /  YAC 

XIX  -  ZZET  •  YAC 

ZETAX-  -  ZXI  •  YAC 

XIZ  -  -XIFT  *  YAC 

ZETAZ-  XX I  •  YAC 


AX  -  AXI  •  XIX  +  AZET  •  ZETAX 

OZ  -  UXI  •  XIZ  4  UZET  *  7ETAZ 

n  •  VXI  •  XIZ  •*  VZET  •  ZETAZ 

AZ  -  AXI  •  XIZ  *  AZET  •  2ETAZ 

TXX  -  -(-4.  •  VX  4  2.  *  VZJ  •  CNM  3. 

TXY  -  CNF  •  (UZ  4  VX) 

•nry  -  -ON*  /  3  •  (-4 .  *  Vz  4  2 .  •  CX) 

K4  •  (<V(I  R)4U(I,«-1))*TXX4<V(I,X)4V(I,F-1))*TXY)*0.5 
1  4  CN.M  /  PP/ (GAMMA  -  l  )  *  AX 

54  -  ((l<I,f)4U<I,X-M)*TXY4<V(I,I)4V(I,K-l))4TYY)*0.5 

1  CAM  /  PR  /  ( GAMMA  -  1.)  *  AZ 

R4  -  0 
S4  -  0 
RH1(K)  -  0 

RH2(F)  -  (ZETAX  •  TXX  4  ZETAZ  •  TXY)  /  YAC 
RH3 ( X )  •  (ZETAX  •  TXY  4  ZETAI  *  TYY )  /  YAC 
60  RH4 ( X )  •  (ZETAX  •  R4  4  ZETAZ  •  S4 )  /  YAC 
DO  TO  F  -  2  ,  KMAXM1 

DOKI.I)  -  DQl(I.F)  -  RBI  (  F  +  l )  -  RB1(F) 

D02(I,R)  -  DQ2 ( I  , X )  ♦  RH2 ( F-i )  -  RR  2 ( F  ) 

DO';I,F)  *  DOJ(I.F)  •»  RBJ(X-l)  RHJ(F) 

DQ4(i,F)  *  DQ4 ( I , X )  *  RB4(K-1)  -  RB4 ( F  ) 

70  CONTINUE 


SOURCE  PROCRAM 

nse2d.f 

time  1 1:46:04  am 

LINE# 


1757 

1758  I  C* 
11759  f  c« 

1760  j  C* 
T76I  !  c- 
1762  '  c« 

1763 

1764 

1765 
"1766 

"1767 
1768  c 
"1769  C 

1770  c 

1771  c 

1772 
"1773 
"1774 
"1775 
1776 
"1777 

1778 

1779 

1780 

1781 
T782 

1783 

1 784 

1785 

1786 


SOURCE  TEXT 


SUBROUTINE  LOW) ( C  L ,  CDP ,  CDF ,  CM  ALFAS  , XREF) 


SUBROUTINE  LOAD 


PARAMETER  ( IX-UO  ,kx-60 1 
C0*«0N'GRID1/X<IX,XX1 ,Y(I1,IX) 

C0F040N/SEINCF/CF ( IX ) 

COHHON/DGRI D/DT , I MAX , EKAX , ITEL , ITEU 
COKMON/SUR  F/PSUR  < IX ) 

TNIS  SUBROUTINE  COMPUTES  TBS  INVISC1C  COOTRIBUTIWS 
TO  LOADS  ON  THE  AIRFOIL  SURFACE 

CL  -  0. 

CD  -  0. 

CDF  -  0. 

CM  -  0. 

DO  300  I  -  ITEL  ,  ITEU  -  1 
DX  -  1(1-11,1)  -  X(I  ,1) 

300  CDF  -  CDF  *  (  CF(I)  *  CF(1U)  )  •  0.5  «  DX 
DO  400  I  «  ITEL  ,  ITEU  -  1 
XL  -  .5  •  (X(l ,1)41(141,1)) 

XL  -  .5  •  (Y(I,1)+Y(I+1,1)) 

DX  -  1(1*1, 1)  -  1(1,1) 

DY  -  Y(I*1,1)  -  Y ( I ,1 ) 

CPA  -  PSUR(IAl)  •  .5  ♦  PSUR(I)  *  .5 
DCL  -  CPA  •  (“DX ) 

DCD  -  CPA  •  DY 


CD  *  SIN(ALFAS) 


1787 

CL  - 

CL  ♦ 

DCL 

1788 

CD  - 

CD  ■» 

DCD 

1789 

400  CM  - 

CM  + 

DCD 

1790 

C 

1791 

DCL 

-  CL 

COS 

1792 

CDP 

-  CL 

SIN 

1793 

CL  « 

DCL 

1  794 

RETURN 

1795 

END 

SOURCE  PROCRAM 

Q3  S3 1 

nse2d.f 

1 1:46.-04  am  | 

PACE* 


24 


SOURCE  TEXT 


SUBROUTINE  TABINT(XP ,YP .XSING , YSING ,N } 


SUBROUTINE  TABINT 


PARAMETER  (IX-180 , Xx-60) 

DIMENSION  XP(IX),YPUX)/SO(IX)  ,AO(IX) 

C  -  XP<1)  “  XSING 

V  -  YP(1)  "  YSING 
U  -  X. 

V  -  0. 

ANGL  -  8.  *  ATAN(l-) 

DO  1  I  -  1,N 

111  -  IP (I)  -  XSING 

Til  *  YP ( I J  -  YSING 

ANGL  *  ANGL  4  ATAN2 ( <U*Y11-V*X11) , (U'XlUV Y1 I ) ) 


2) 


5) 

-5) 


R  -  SQRT(X11**2  4  Yll 
0  -  XU 
V  -  Yll 
R  -  SQRT(R) 

A0( I )  -  R  •  COS (ANGL 
S0(I)  -  R  •  SIN( ANGL 
DX  -(A0(N^-A0(l))/96. 

AOST  *  AO ( 1 ) 

DO  1  I  •  1  r  »? 

XX  -  FLOAT(I-l)  *  DX  4  AOST 

CALL  TAINT ( AO  ,S0  ,  XX  ,  YY  *N  ,  3 ,NER ,  MON ) 

XP(I)  -  XX  -  XX  -  YY  *  YY  ■»  XSING 

YP  ( I )  -  2.  •  XX  •  YY  4  YSING 

RETURN 

END 


I 


SOURCE  PROCRAM 

nse2d.f 


DATE  1/20/90  PAC£# 
^nme  1  U4&04  am  I  25 


SOURCE  TEXT 


SUBROUTINE  T8INT(XTAB,rTAB,X,rX,N,I,NER,MON) 


1  874  C* 

1875  C* 

1 876 
T877 

1878  C 

1879  c 

1880  c 

1881  c 

1882 
1883 

”1884 

1885 

1886 
TS87 
1888 
1889 

”1890 
”189 1 
T892 
'  1893 


M BROUTINE  TAINT 


PARAMETER  (IX*l«0,kx*',0) 

DIMENSION  ITAB<1) ,FTAB(1) ,T(10) ,C(10) 

mass  -  AXES  SUBROUTINE  POX  POLTNOMI M.  INTERPOLATION 
Or  A  TABULATED  FUNCTION 

ir(N-X)  1,1,2 

1  NER  -  2 
RETURN 

2  ir(X-»)  3,3,1 

3  ir(MON)  4.4,5 

5  ir<NON-2>  6,7,4 

4  J  -  0 

NM1  *  N  -  1 
DO  •  I  *  1  ,  NK1 
ir(XTAB(I)  -  XTABt I+1 ) )  »,li,10 
11  NER  -  3 


1894 

9 

J  -  J-l 

1 89  S 

GO  TO  8 

1896 

10 

J  ■  J+l 

1897 

8 

continue 

1898 

MON  -  1 

1899 

IF<J)  12  ,  6  ,  6 

1900 

12 

MON  -  2 

1901 

7 

DO  1J  I  •  1  ,  N 

1902 

IT(X  -  XTAB(I))  14,14, 

1903 

14 

J  -  I 

(904 

GO  TO  18 

190$ 

13 

CONTINUE 

1906 

GO  TO  15 

1907 

6 

DO  16  X  -  1  ,  N 

”"1908 

IF(X-XTAB(I) )  16.17,17 

1909 

17 

J  -  I 

1910 

GO  TO  18 

191  1 

16 

CONTINUE 

1912 

15 

J  ■  N 

1913 

18 

J  -  J  -  (K  +  l)  /  2 

1914 

I  F( J )  19 ,  a  9 , 20 

191$ 

19 

J  *  1 

17'  J 

20 

M  -  J  +  X 

1917 

IF(K  -  N)  21,21,22 

1918 

22 

J  -  J  -  1 

1919 

GO  TO  20 

1920 

21 

JTPl  -  X  +  1 

1921 

JSAVE  -  J 

1922 

26 

DO  23  1  -  1,  KP1 

1923 

C ( 1 )  -  X  -  XTAB(J) 

1924 

T(L)  -  FT^BfJJ 

1925 

t  23 

J  -  J  +  l 

1926 

, 

DO  24  J  »  1,X 

1927 

! 

I  -  J  +  l 

1928 

25 

T(  I  1  -  (C(J)*T(I)-C<I) 

1929 

I  -  1  +  1 

1930 

IF(I-XPl)  25,25,24 

1931 

2 4 

CONTINUE 

1932 

FX  -  T(KP1) 

1933 

NEP  -  1 

1934 

RETURN 

193$ 

END 

SOURCE  PROCRAM 


date  1/20/90  I  *ace# 
t>me  11:464m  am  i  ‘ 


SOURCE  PROGRAM 

nse2d.f 


DATE  1/20/90 
~™«  11:4604  am  2S 


SUBROUTINE  CLVStt(D$ ) 


SUBROUTINE  CLUSTR 


PARAMETER  (IX-180,k.x-60) 

COmOH/GR  I  Dl/X  (ll,EX)/E(IX/RX) 

COM80W/DCRI D/DT ,  IMAX ,  RMAX ,  I  TEL ,  I TEU 
DIMENSION  S { 6 0 ) , XP (60), YP ( 6 0 ) , R ( 6 0 ) 

T1IS  SUBROUTINE  CLUSTERS  A  GIVEN  X,Z  GRID  SUCfi  THAT  TEE  FJRST  POINT  IS  AT 

DO  100  I  -  1  ,  I MAI 
S(l)  -  0. 

XP(1)  -  X(I,1) 

TP(1)  -  Z(l,l) 

DO  10  I  -  2  ,  RMAX 
XP(R)  -  X( 1 , K  j 
TP(R)  -  Z ( I ,  R ) 

10  S (X)  *  SQRT((XP(R)-XP(R-1))**2+(YP(R)-YP(R-1))**2) 

1-*S(R-1) 

SUKDX  -  S(KMAX) 

CALL  STFTCH ( SUMDX  , DS , El , RMAX , FACTOR ) 

WRITE  (8,200)  1, FACTOR 
-  0. 

D*  -  DS 

DO  20  X  -  2  ,  RMAX 
R(R)  -  R(R-l)  *  DR 
DF  »  DR  •  k ACTOR 
20  CONTINUE 

FLAST  -  1.  /  H(RMAX) 

DO  30  R  -  2  ,  RMAX 

R1  »  t(R)  *  FLAST  *  S(KMAX) 


DO  30  R  -  2  ,  RMAX 

R1  »  R<  R )  •  FLAST  *  S(KMAX) 

CALL  TAINT < S , XP ,R1 ,XP1 , RMAX , 3 ,NER ,MON ) 

X(I ,  R)  -  XP1 

CALL  TAINT ( S , YP,Rl , YP1 , RMAX , 3 , HER, MON ) 

Z(I,R)  -  YP1 
30  CONTINUE 
100  CONTINUE 

WRITE  (6,115) 

DO  110  1*1,  IMAX 
DX  *  X( I , 2 )  -  X(I,1) 

DY  *  Z(I ,2)  -  Z(1 ,1) 

DN  -  SQR  T  ( DX  *  DX-*  DY  *  DY ) 

WFITE( 6,120)  I  ,  DX  ,  DY  ,  DN 
110  CONTINUE 
RETURN 

115  FORMAT( 5X , 6HNORMAL, IX , 8HDISTANCE , 3H  AT , 4H  THE , 5R  WALL,/ 
1, 5H  I , 8X , 2HDX , 8X , 2HDY , SX, 2BDN ,// ) 

120  FORMAT (15, 3E10 . 5) 

200  FORMAT ( 1 5 , FI 0 . 5 ) 

END 


Wm **  fsftjtsi  WwllS/WW m|m' :W^Mi*wWW  NN *«\**jWl  C 


139 

1 40 

1 4  r 

142' 
1431 
I  44 

145 

146 

147 
148' 

149 

150 

151 
152' 
153 

I  54' 
155 


156 

is  r 

158, 

1591 

l«T 


161 

162 

163 

164 

165 

166 

1 67 

168 

169 


SUBROUTINE  STRTCB<SUMDX ,DXI , Tl ,N1 ,F ) 

C* 

C*  SUBROUTINE  STRTCH 

C' 

c*. *«%«•*•«*«««••«*»**«•***•*«««•<•* v************ 

PARAMETER  <IX-180,ltx-60) 

C 

C  ms  SUBROUTINE  DETERMINES  A  GEOMETRIC 

C  PROGRESSION  OF  GRID  SPACING  BETNEEN  1  AND  hi  SUB  THAT 
C  SUNUDX}  EQUALS  SUJOI.  TIE  RATIO  BETWEEN  SUCCESSIVE 
C  BRACINGS  IS  R, 

N  -  VI  -  1 
R  -  1.5 
El  -  l.E-04 
E2  -  l.E-04 
DO  10  L  -  1,  50 

P-  (R-l )  *  SUMDX  -  DX1*(R*«N-1) 

FP  -  SUMDX  -  DX1  •  FLOAT(N)  *  R  **  (N-l) 

RITER  *  P  -  F/  FP 

C  IF( l.E-02. IT. RITES. AND. RITER. LT. 1. >  RITER  -  1. 

C  TE(1. .LI.FlTfiR.AND.RITER.LT.  100.  )  RITER-. 01 

IF( ABS (R-RITER J . LT.  R*E1)  GO  TO  1 
R  -  RITER 
10  CONTINUE 
R  -  1.0001 

C  DU  -  OTTOT/FLOAT(N1-1) 

RETURN 
1  R-  RITER 
RETURN 
END 


f 


i 


I 


SOURCE  PROCRAM 

nse2d.f 


DATE  1/20/90  fMl* 

11:46904  am  I  ^0 


SOURCE  TEXT 


SOMOOTXN E  EDDY  ( COT L) 


1173  i  c* 

l£2f  c* 


SUWOOTIW  EDDY 


t  I  BH 

IT*  i 

~fi*2 

~2T84'j 
2f«s1  c 


2(98 
2199  c 
'2200 
32  01 
2202 

2203 

2204 

2205 

2206 

2207  ! 

2208 


PARAMETER  ( IX-180 , »UC«60) 

COI«ON/FLON/Ql(IX,n>,Q2(IX,ia)  ,Q3(IX,KX)  ,Q4(IX,KX) 
COHMON/MUTUR/CHU ( IX , XX ) 

OOMION/SX  INCF/Cr(  ix  ) 

COMMON /DGP.  I D/DT , I MAX , XMAX , I TEL ,  ITEO 

COMMON /PAP/GAMMA , REYREF , ALFA, ALTAI /REDFRE , AMINF , ALTAI 
COMION/GRI Dl/X ( IX ,  EX  ) , Z ( IX , XX ) 

COMMON/SPEED/  U(IX,KX) ,V(IX,KX) ,AA(IX,XX) 

DIMENSION  TI N  ( EX ) ,  TOUT  ( KX  )  ,  Y  ( XX  ) ,  WANS  ( I X )  ,  S  ( I X )  ,  UU  (  XX ) ,  UE  (IX) 

CMUPP  -  1000.  '  CMUL 
KEDGE  -  KMAX 

ILE  -  (  ITEL  •»  ITEO  )  /  2 
CHORD  -  X  ( ITEV ,  1 )  -  X(ILE,1) 

DO  100  I  ■  2  ,  IMAX  -  1 

ir<  AB£(X( 1,1)). LT.  <  AB5(X(1LE,1) )  ♦  0.05  ♦  CHORD  )  J  GO  TO  100 
UDIF  -  0. 

UMAX  -  0. 

OMIN  -  9999. 

FMAX  -  0 . IE- 06 
YHAX  -  . IE-06 
PYMAX  -  0. IE-06 
Y<1>  -  0. 

COMPOTE  TAD  AT  THE  HALL 
UET  -  0(1.2)  -  U(X,l) 

VET  -  V( I . 2 )  -  V(I,1) 

XXI  -  X(X~1,1)  -  X(I-1,1) 

2X1  -  2<IU.l)  -  2(1-1. 1) 

XET  -  4.  •  Ml, 2)  -  3.  *  X(l,l)  -  1(1,3} 

ZET  -  4.  «  2(1,2)  -  3.  *  2(1,1)  -  2(1,3) 

XXI  -  .5  •  XXI 

2X1  3  .5  •  ZXI 

XET  -  .5  *  XET 

ZET  -  .5  •  ZET 


7210 

YAC  -  1.  (XXI  *  ZET  -  ZXI  *  XET) 

221 1 

OMEGA  -  (UET  •  XXI  +  VET  *  ZXI  )  •  YAC 

2212 

TKALL  -  A.V1NF  •  OMEGA  /  REYREF 

2213 

CF ( 1 )  -  2  •  TWALL  /  ( AMINF*  *  2 ) 

2214 

FACT  *  SOr,T(Ql(I,l)  •  AB5  (  TWALL)  J  .REYREF/  (  26  -  *  AMINF) 

2215 

DO  10  K  -  2  ,  KEDGE- 1 

2216 

UXI  «  U ( 1 -  1 , K )  -  U ( 1 - 1 , K ) 

2217 

VXI  -  V(I-l.K)  -  V ( 1 - 1 , F ) 

2218 

UET  -  U(I  .  K-*  1 )  -  U(I  ,K-1) 

2219 

VET  «  V(Z  ,  KU)  -  V(I  ,X-1  ) 

2220 

XXI  -  X ( I- 1 , K )  -  X(X-l.K) 

2221 

ZXI  »  Z(l- 1  ,K)  -  Z ( I - 1 , K ) 

2222 

XET  -  X ( I .K41)  -  X(I,K-1) 

2223 

ZET  -  Z(I ,  K  + 1  )  -  2(1, R-l) 

2224 

YAC  -  1.  /  (XXI  *  ZET  -  ZXI  *  XET) 

2225 

OMEGA  -  ABS(UET*XX1-VET«ZXI-UXI*XET-VXI*ZET)  *  YAC 

2226 

UTOT  -  SORT  (  U  ( I  ,  K  )  •  *  2  •*  V<I,K)*«2) 

2227 

UMAX  -  AMAX 1  ( UTOT  ,  UMAX  ) 

2228 

UMI N  •  A'-!  I K 1  (  UTOT  ,  UMI N  ) 

222 9 

Y(K)  -  S©fcT( (X( X , K)-X( I , K~1 ) ) * *2+ ( 2( I ,K)-Z( 1 ,K-1 ) ) * • 

2230 

F  -  Y(K)  -  OMEGA 

2231 

IF( ( Y ( K ) *  FACT) . GT . 20 . )  GO  TO  31 

2232 

IF(I  CT.ITEL.ASD.I .LT.ITEU)  F  «  F  *  (1.  -  EXP(-Y(K)« 

2233 

31 

CONTINUE 

2234 

c 

2235 

I F ( F . GT . FYAX )  THEN 

2236 

FMAX  «  r 

2237 

YKAX  “  Y ( K ) 

2238 

ENDIF 

2239 

FCT  -  Y(X)  •  FACT 

2240 

I F ( FCT . GT  20. )  FCT  -  20. 

2241 

FCT  -  ABS(FCT) 

2242 

EL  -  .4  •  Y ( X )  *  (1.  -  EXP ( -FCT) ) 

2243 

TIN ( X )  -  Q1 ( I , K )  *  EL  •  EL  *  OMEGA 

2244 

TIN(X)  -  ASS ( TIN( X ) ) 

2245 

10 

CONTINUE 

2246 

UDIF  -  ABS ( UMAX 'UM IN ) 

2247 

KSWTCB  -  0 

2248 

FWAKE  -  YMAX  *  FMAX 

2249 

FI  -  0.25  •  YMAX  *  UDIF  **2  /  FMAX 

2250 

i 

IF(F1 . LT. rWAKE)  FWAKE  -  FI 

2251 

i 

DO  20  X  -  2  ,  KEDGE  -  1 

2252 

j 

FKLEB  -  0 

2253 

IF(ABS(Y(K)/YMAX)  LT.l.E+04)  THEN 

2254 

1 

rXLEB  •  1  /  (1.  +  5.5  •  (0.3  «  Y(K)/YMAX)  ••  6) 

2255 

| 

END  IF 

2256 

TOUT ( K )  -  .0168  *  1.6  •  Q1(X,K)  «  FKAKE  •  FKLEB 

2257 

TOUT <  K)  -  ABS ( TOUT ( X ) ) 

2258 

IF(KSWTCB  NE.O)  GO  TO  20 

2259 

I 

I F ( TI N ( X ) . GT . TOUT ( X ) )  ISNTCB  -1-1 

2260 

20 

CONTINUE 

2261 

DO  30  X  •  2  ,  KEDGE  -  1 

2262 

IF(  K  .  LE  KSWTCB  )  THEN 

2263 

CMU ( I , K )  -  ABS(TIN(K) ) 

2264 

ELSE 

2265 

CMU ( I , K )  -  ABS(TOUT(K)  ) 

2266 

END  IF 

2267 

30 

CONTINUE 

2268 

C 

2269 

c 

molt  TIE  VALUE  OF  EDDY  VISCOSITY  TO  AN  UPPER  LIMIT 

2270 

c 

m\ 

c 

i  -  l 

2272  !  C  734  R  -  X  4  1 

2273  i  C  irtX.CT. HEDGE)  GO  TO  736 

2274  j  c  735  IF(CMV(J,I)  LE. CMUPP}  GO  TO  734 

2275  c  cmu ( i , x )  -  cmupp 

2276  t  c  i  -  x  ■*  l 

2 277  !  C  IF( K. LE  XEDGE )  GO  TO  735 

2276  !  C  736  continue 

2279  100  CONTINUE 

2260  '  RETURN 

2261  ‘  end 


SOURCE  PROGRAM 

nse2d.f 


DATE 


1/20/90 


11:4&04am 


FACE  • 


31 


UNEf 


SOURCE  TEXT 


2  282 
p 

2284 

fees' 

2286 

J287 

12288 

^2289 

2290 

2291 

2292 

2293 

2294 
2  29S 

2296 

2297 

2298 

2299 
.2306 
£2301 
£2302 

2303 

2304 

2305 

2306 

2307 

2308 
12309 
"2310 
'231 1 
2  3 1 2 

2313 

2314 

2315 

2316 
231  7 

2318 

2319 

2320 

2321 

2322 

2323 

2324 

2325 

2326 

2327 

2328 

2329 

2330 

2331 

2332 

2333 

2334 

2335 

2336 

2337 

2338 

2339 

2340 

2341 
234  2 

2343 

2344 

2345 

2346 

2347 

2348 

2349 

2350 


SUBROUTINE  RTS I 


SUBROUTINE  HE  SI 


PARAMETER  (IX-180 ,fcX-«0) 

COMHON/r I X /OMEGA , HOOT 

CQM*DN/PERTR/DQ1(II ,  EX ) , DQ2 ( IX ,  XX  )  ,  DQ3 ( IX  ,  XX )  ,DQ4(1X.XX) 

COW40N/GR I  Dl/X  ( IX ,  EX ) ,  Z  <  IX ,  IX ) 

COM40N/DGR I D/DT ,  IMAX ,  KMAX , I TEL , I TEU 
CO*#«ON/FLOW/Q1  { I X ,  EX ) ,  Q2  ( IX ,  XX } ,  Q3  ( IX ,  IX ) ,  Q4  (  I X ,  XX ) 

COMION/SPLED/  U  ( IX ,  XX )  ,V(  IX  ,XX )  ,  AA(  IX ,  XX  ) 

COHMON/PAR/G AMMA ,  REYREF ,  ALFA ,  ALFA1 ,  REDFRE ,  AMI  NT,  AX  FA  2 
C0M40N/C0KRHS/  RHS (IX, 4) 

XTAU<2 /X)  -  OMEGA  •  2(1, X) 

YTAC ( I ,  X )  -  -  OMEGA  •  X(I,X)-  ROOT 

THIS  SUBROUTINE  COMPUTES  THE  RESIDUAL  ON  THE  RIGHT  HAND 
SIDE  ARISING  FROM  THE  EULER'  PART  Of  THE  GOVERNING  EQUATIONS 

FLUX  TERMS  NI THIN  THE  IS-  fiCRIVATIFE 
DO  100  X  -  2  ,  XMAX  -  1 
DO  10  I  -  1  ,  IMAX 

UCON  -  U(I,X)  •  (2(I,X+1)-Z(I ,X-1) ) 

1  -  V(I,X)  *  (X(I , X+I)-X(I ,X-1) ) 

UCON  «  0.25  *  DT  •  UCON 

XI T  -  -  XTAU ( I  ,  X )  *(2(I ,X+1)-Z(I ,  X-l ) ) 

1  4  YTAU  ( I  ,  X )  «  (X(I,K*1)  -  X  ( 1 ,  X- 1 )  ) 

XIT  -  XIT  •  DT  *  0.25 

UCON  -  UCON  -p  XIT 

RHS ( X  ,  1 )  -  UCON  •  01<I,X) 

P  -  (GAMKA-I.)  *  ( Q4 ( I , X )  -  . 5-Q1 < I , X ) • < U < I , X ) • *2  4  V(I,X)**2)  ) 
RHS (1,2)  «  Q2(I,X)  *  UCON  *  P  *  DT  •  0.25  *  <2(3, X+l)  -  Z(I,X-1)> 

RHS( I , 3 )  «  Q3(I,X)  *  UCON  -  P  •  DT  *  0.25  •  ( X ( 1 , X+l )-X ( I , X-l ) ) 

RH$(I,4)  -  UCON  *  (Q4 (1 , X)*P)  -  XIT  •  P 

10  CONTINUE 

DO  11  I  -  2  ,  IMAX  -  1 

DQ1 ( I , X )  -  DQ1 { X  ,X)  -  RHS( 1*1 ,1)  ♦  RHS(I-1,1) 

DQ2 ( I , X )  -  DQ2 ( X , X )  -  RHS (1*1,2)  4  RHS(I-1,2) 

DQ3 ( I , X J  -  DQ3(I,X)  -  RHS (1+1,3)  4  RBS(I-1,3) 

11  DQ4 ( I , X )  -  DQ4 ( I , X )  -  RHS(l4l,4J  4  RHS(I-1,4> 

100  CONTINUE 


FLUX  TERMS  WITBIN  THE  ETA-  DERIVATIVE 

DO  200  I  -  2  ,  IMAX  -  1 

DO  20  X  -  1  ,  XMAX 

VCON  -  U(J,K)  «  <Z<Z-1,X>-2(X+1,X)) 

1  +V(I,X)  *  (X(I+1,X)-X(I-1,X)| 

VCON  »  VCON  *  0.25  •  DT 

ETAT  -  -XTAU(I,X)  *  ( Z< I- I , K >-Z ( 1+1 . X ) )  -  YTAU(I,K)* 

1  <X(X*1,X)-X(I-1,X)  ) 

ETAT  *  ETAT  *  0.25  •  DT 

VCON  -  VCON  +  ETAT 

RHS ( X , 1 )  •  VCON  •  Q1(I,X) 

P  -  (GAMMA-1.)  •  ( Q4 ( I , X )  -  0.5  •  Ql( I . X) • (U( I , X) *«2  +  V(1,X)*«2)) 
RHS (  X  ,  2 )  *  VCON  •  Q2  ( I  ,  X  )  +P  •  DT  •  .25  «  < 2 < 1-1 , X )-Z < I +1 , X ) ) 

RHS ( X , 3 )  -  VCON  »  03(1, X)  +  P  *  DT  *  . 25  *  (X(Z*1,X)  -  X(I-1.X)J 
RHS ( X , 4 )  -  VCON  •  (Q4(I,K)*P)  -  ETAT  •  P 
20  CONTINUE 

DO  21  X  -  2  ,  XMAX  -  1 
DQl(I.K)  *  DQ1 ( I , X )  -  RBS(X+1,1) 

RHS ( X*1 ,2) 


DQ2(I,X) 


DQ2 ( I , X ) 

DQ3 ( 1 , X )  «=  DQ3 ( I , X )  -  RBS(X+1,3) 
21  DQ4 ( I , X )  «  D04(I,K)  -  RHS ( K*1 , 4 ) 
200  CONTINUE 
300  FORMAT (21 6 ,4E14 . 4) 

RETURN 

END 


RHS ( X- 1 , 1  ) 
RHS ( X- 1 , 2 ) 
RHS ( K-l , 3 ) 
RHS (X- 1,4) 


i 


IRCE  PROCRAM 


1/20/90  I  PACE# 


t»«e  11:46.-04  am 


SUBROUTINE  R0TGRID(  IMAX ,  KMAX ,  DALFA) 


SUBROUTINE  ROTGID 


PARAMETER  <IX-180,Rx-60) 

ROTATE  GRID  IN  THE  CLOCKWISE  DIRECTION  BY  AN  AMOUNT  DALTA 
COHMON/GR1D1/X{IX,KX) ,Z{IX,KX) 

CA  -  COSiDALFA) 

SA  -  -  SIN(DALFA) 

DO  20  K  -  1  ,  KMAX 
DO  20  I  -  1  ,  IMAX 
XI  -  I(1,K) 

Z1  -  Z(1 ,X) 

X{1 , K)  »  XI  •  CA  -  El  *  SA 
20  Z(I , K)  -  Z1  *  CA  4  XI  •  SA 
RETURN 
END 


SOURCE  PROGRAM 
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DAT! 


1/20/90 


T,ME  1 1:46XM  am 


PACE* 


33 


UNE  # 


SOURCE  TEXT 


2170 

2371 

2372 
7373 
2374 
1375 
ff  76 
2377 

7378 

7379 

7380 

2381 

2382 
7383 
2384 
2365 

"7366 

7367 

7368 

2389 

2390 

2391 

2392 

2393 
7394 

2395 

2396 

2397 

2398 

2399 

2400 

2401 

2402 

2403 

2404 

2405 

2406 

2407 

2408 

2409 

2410 
241  i 


SUBROUTINE  CPPLOT(Il,l2  , FMACB,X,Y) 


PARAMETER  ( IX-1B0, Ax-60) 

THIS  SUBROUTINE  PLOTS  CP  AT  EQUAL  INTERVALS  IN  TEE  MAPPED  PLANE 
COMMON/ SURF/PS UR (IX) 

DIMENSION  RODE  (4  )  , LINE <  90)  ,  X(  IX  )  ,  If  (IX) 

DATA  KODE/IH  , IB* , 1HI , 1H •/ 

NRXTE  (4,2) 

2  FORMAT  (5  OB  OP  LOT  OF  CP  AT  EQUAL  INTERVALS  IN  THE  MAPPED  PLANE/ 

1  10H0  X  ,  10H  X/C  ,  1  OH  CPL  ,10B  CPI  > 

CPO  -  (1.  4  .2  *  FMACH  *  *2 )  *•  3.5  -  1. 

CPO  -  CPO  /  (  .7  •  PMACH  * *2 ) 

10- 30.  •  CPO  *  4.5 
I  MIN  -  (I2~U)/2  4  II 
ILOH  -  2  *  I  MIN 
CHIHX(II)  -  X(IMIN) 

DO  12  1  -  1  ,  90 

12  LINE(l)  -  KODE(l) 

DO  14  I  -  I MIN  ,  12 

1-30.  *  (CPO  -  PSUR(If)  *  4.5 

11- 30.  *  (CPO  -PSUR(ILOW-I))  4*5 
IF(X.LT.l)  X  -  1 

IFCK1.LT. 1)  II  -  1 
1 F( X . GT . 90 )  1-90 
IF( XI  .GT.  90)  XI  -  90 
LINE(KO)  -  KODE ( 3 ) 

LINE ( X )  -  XODE ( 2  J 
LINE(Xl)  -  XODE ( 4 ) 

XOC  -  <X(I)  -  X i TMIN ) )  /  CBD 

NRITE  (6,610)  X ( I ) ,XOC , PSUF ( I IOW-I ) , PSUR ( I ) , LINE 
LINE ( XI )  -  XODE ( 1 ) 

34  LINE ( X )  -  XODE ( 1 ) 

RETURN 

600  FORMAT ( 4 FI 0 . 4 ) 

610  FORMAT ( 4F1 C . 4 , 90A1 ) 

END 


SOURCE  PROGRAM 

nse2d.f 


-*?*-*m*M 


DATE 


1/20/90 


™£  1 1:46*>4  am 


PACE* 


34 


LINE  * 


SOURCE  TEXT 


2412 
241  3 
>414 


SUBROUTINE  OUTRVC  ( REPEAL ) 


24 1 C 
>416 

241 7 

2418 

2419 

2420 

2421 

2422  j 

2423  i 

2424  ! 

2425 

2426 
>427 
2426 

2429 

2430 

2431 

2432 

2433 

2434 

2435 

2436 

2437 

2438 

2439 

2440 

2441  ; 

2442 

2443  { 

2444 

2445  ! 

2446  ' 

2447 ; 

2448 

2449 

2450 


SUBROUTINE  OUT* VC 


PARAMETER  iIX-lSO.Ax-SO) 

COIWOK/PLOT/TITLE(  10 )  ,NSTPT ,  RESD(  3000)  ,RES,CLB<  3000)  ,CDPB( 3000  ) 
COM40N/DGF I D/DT ,  IMAX ,  MAX ,  I  TEL ,  I  TEL’ 

CO»40N/rLON/Ql<IX,KX),Q2(IX,  KX )  .03  (IX,  IX)  ,Q4(IX,KX) 
C0HM0N/PAR/6AMMA , REYREF , ALFA , ALFA1 , REDFRE , AMINF , ALFAI 
COKMON/SURF/PSCR<IX) 

COMMON/SKI NCF/Cr<  IX) 

COMMON/MU  TUR/CMl)  ( IX ,  XX  ) 

COMMON/TITL/I  TITLE 
CHARACTER  ITITLE'80 
OINENF  ON  EY<  IX ,  KX  ) 

DOM  -  0. 

CODE  •  0. 

IRES  -  IFIX(RES) 

AIT  AD  -  ALFA  •  4  5.  /  ATAN(l-) 

CKOL  -  AMINF  /  REYREF 
DO  10  K  ■  1  ,  KMAX 
DO  10  I  •  1  ,  IMAX 

nr<i,K)  -  o. 

10  CM0(I,R)  •  CMU(I.I)  /  CKL'L 
WRITE ( 4  )  1TITLE 

WRITE (4 )  IMAX, KMAX , ITEL, ITEU , AMINF, ALFAD, REREAD, DDM,NSTPT, GAMMA, 
♦  CODE , RES . DUM 
WRITE ( 4 )  Q1 , QT , Q3 , Q4 
WRITE ( 4 )  RESD 
WRITE ( 4 )  PSUR 
WRITE  (  4  )  CF 
WRITE  (  4  )  CLB.CDPH 
WRITE ( 4 )  CMC 
WRITE {4 )  FT 
WRITE( 4 )  TK 
WRI TE ( 4 )  TE 

return 


I 


SOURCE  PROGRAM 

date  1/20/90 

nse2d.f 

^  1 1:46.04  am 

5£.  M*.  ...=s£2fe''&s££i: 


SOURCE  TE>rr 


2452 

SUBROUTINE. 

OUTOT(NOUT, NSTPP, TIME 

REREAL) 

2453 

PARAMETER 

(IX-l«0,*x-60) 

2454 

COMMON/CRI Dl/X { IX , KX ) , Z ( IX  ,  KX } 

2455 

COMMON/PAR/GAMMA ,  REYREF ,  ALFA ,  ALFA1 ,  REDFTtE  ,  AMI NF ,  ALFAI 

2456 

COMMON/DGRID/DT ,  IMAX ,  KMAX  ,  I  TEL,  ITEU 

2457 

COMMON /FLOW/01  ( IX , KX )  ,Q2  ( IX  ,  KX  )  ,Q3  ( IX  ,  EX )  ,Q4(IX,KX) 

2458 

ITN  •  NSTPP 

2459 

PI  »  4 . *  AT AN ( 1 . ) 

2460 

ALFAD- ALFA 

* (180. /PI ) 

2461 

IF( 

ITN  .EQ.  1 *NOUT  )  THEN 

2462 

REWIND  31 

2463 

WRITE  <  31 ) 

IMAX  ,  KMAX 

2464 

WRITE  (31) 

(  (  1(1,1),  1-1 , XMAX 

,  1-1, KMAX  ) 

2465 

WRITE  (31) 

(  (  1(1,1),  I -1, IMAX 

,  1-1, KMAX  ) 

2466 

WRITE  (31) 

IMAX  ,  KMAX 

2467 

WRITE  (31) 

AMINF,  ALFAD,  REREAL, 

TIKE 

2468 

WRITE  (31) 

<  (  01(1, X),  1*1, IMAX 

),  1-1, KMAX  ) 

2469 

WRITE  (31) 

<  (  02(1,1),  1-1, XMAX 

),  1-1, KMAX  ) 

2470 

WRITE  (31) 

(  (  03(1, X),  1-1, XMAX 

),  1-1, KMAX  ) 

2471 

WRITE  (31) 

(  (  04(2,1),  2-2 , XMAX 

),  K  *1 ,  KMAX  ) 

2472 

REWIND  31 

2473 

END  IF 

2474 

IF( 

ITN  .EQ.  2«NOUT  )  THEN 

2475 

REWIND  32 

2476 

WRITE  (32) 

IMAX  ,  KMAX 

2477 

WRITE  (32) 

(  (  1(1,1),  1*1, XMAX 

,  1-1, KMAX  ) 

2478 

WRITE  (32) 

(  {  Z(I,K),  1-1, IMAX 

,  1*1, KMAX  ) 

2479 

WRITE  (32) 

IMAX  ,  KMAX 

2480 

WRITE  (32) 

AMINF,  ALFAD,  REREAL, 

TIME 

2481 

WRITE  (32) 

(  (  Q1 ( 1 , K } ,  1-1, IMAX 

),  1-1, KMAX  ) 

2482 

WRITE  (32) 

(  (  02(2,1),  2-2, IMAX 

),  1-1, KMAX  ) 

2483 

WRITE  (32) 

(  (  Q3 ( I , 1) ,  1-1, IMAX 

),  1-1, KMAX  ) 

2484 

WRITE  (32) 

(  (  04(1,1),  1*1 , IMAX 

),  1-1 , KMAX  ) 

2485 

REWIND  32 

2486 

END  IF 

2487 

3F( 

ITN  .EQ  3 *NOUT  )  THEN 

2488 

REWINE  33 

2489 

WRITE  (33) 

IMAX  ,  KMAX 

2490 

WRITE  (33) 

(  <  1(1,1).  1-1 , IMAX 

,  K-l , KMAX  ) 

2491 

WRITE  (33) 

(  (  2(1,1),  X-1,XKAX 

,  K-l , KMAX  ) 

2492 

WRITE  (33) 

IMAX  ,  KMAX 

2493 

WRITE  (33) 

AMINF,  ALFAD,  REREAL, 

TIME 

2494 

iTITE  (33) 

(  <  Ql(I.K),  1-1,1  MAX 

),  K-l, KMAX  ) 

2495 

WRITE  (33) 

(  (  Q2(I.l),  1-1 , IK*  X 

) ,  K- 1 , KMAX  ) 

2496 

WRITE  (33) 

(  (  03(1, F).  I-l.IM.Ol 

),  K-l, KMAX  ) 

2497 

WRITE  (33) 

(  (  04(1,1),  1-1, IMAX 

),  K-l, KMAX  ) 

2498 

REWIND  33 

2499 

END  IF 

2500 

IF< 

ITN  .EQ-  4 *NOL’T  )  THEN 

2501 

REWIND  34 

2502 

WRITE  (34) 

IMAX  ,  KMAX 

2503 

WRITE  (34) 

(  (  1(1,1),  1-1, XMAX 

,  K* 1 , KMAX  ) 

2504 

WRITE  (34) 

(  (  2(1,1),  I« 1 , IMAX 

,  K-l, KMAX  ) 

2505 

WRITE  (34) 

IMAX  ,  KMAX 

2506 

WRITE  (34) 

AMINF,  ALFAD,  REREAL, 

TIME 

2507 

WRITE  (34) 

(  (  OKU),  1*1,1  MAX 

),  K-l, KMAX  ) 

2506 

WRITE  (34) 

(  {  Q2(I  ,F.)  ,  1*1  ,  IMAX 

),  K-l, KMAX  ) 

2509 

WRITE  (34) 

(  (  03(1,1),  I *1,1 MAX 

),  K-l, KMAX  ) 

2510 

WRITE  (34) 

(  (  04(1 ,K) ,  1=1, IMAX 

),  K-l, KMAX  ) 

25  H 

REWIND  34 

2512 

END  IF 

2513 

I  F  ( 

ITN  .EQ.  5*NOUT  )  THEN 

2514 

REWIND  35 

2515 

WRITE  (35) 

IMAX  ,  KMAX 

2516 

WRITE  (35) 

(  (  X(I ,K) ,  1*1 , IMAX 

,  K-l, KMAX  ) 

2517 

WRITE  (35) 

(  (  2(1 ,K) ,  1-1, IMAX 

,  K-l, KMAX  ) 

2518 

WRITE  (35) 

IMAX  ,  KMAX 

2519 

WRITE  (35) 

AMINF,  ALFAD,  REREAL, 

TINJ 

2520 

WRITE  (35) 

(  (  01(1,1).  1*1 , IMAX 

),  K-l, KMAX  ) 

2521 

WRITE  (35) 

(  (  Q2< I ,K) ,  1*1 .IMAX 

),  K-l.  KMAX  ) 

2522 

WRITE  (35) 

(  (  03 ( 2 , K ) ,  1*1,2 MAX 

),  K-l, KMAX  ) 

2523 

WRITE  (35) 

(  (  04(1,1).  1*1, IMAX 

),  K-l, KMAX  ) 

2524 

REWIND  35 

2525 

END  IF 

2526 

IF( 

ITN  .EQ  6 *N0L’T  )  THEN 

2527 

REWIND  36 

2528 

WRITE  (36) 

IMAX  ,  KMAX 

2529 

WRITE  (36) 

(  (  X(I.K),  1*1,1 MAX 

,  K-l, KMAX  ) 

2530 

WRITE  (36) 

(  (  Z< I , K ) ,  1*1 , IMAX 

),  K-l, KMAX  ) 

2531 

WRITE  (36) 

IMAX  ,  KMAX 

2532 

WRITE  (36> 

AMINF,  ALFAD,  REREAL, 

TIME 

2533 

WRITE  (36) 

(  (OKI  ,  K  )  ,  1*1 ,  IMAX 

),  K-l, KMAX  ) 

2534 

WRITE  (36) 

(  (  Q2< I , K) ,  1-1 , IMAX 

),  K*  1  ,  KMAX  ) 

2535 

WRITE  (36) 

(  <  03(1 ,K) ,  1*1,1  MAX 

),  K-l, KMAX  ) 

2536 

WRITE  (36) 

(  (  Q4(I ,K) ,  1-1, IMAX 

),  K-l, KMAX  ) 

2537 

REWIND  36 

2538 

END  IF 

2539 

IF< 

ITN  .EQ.  7 *NOL'T  )  THEN 

2  540 

REWIND  37 

2541 

WRITE  (37) 

IMAX  ,  KMAX 

2542 

WRITE  (37) 

(  (  X(I  K),  1*1 , IMAX 

).  K-l  ,  KMAX  ) 

2543 

WRITE  (37) 

(  (  2(1,1).  I-l.IMAX 

),  K-l, KMAX  ) 

2544 

WRITE  (37) 

IMAX  ,  IMAX 

2545 

WRITE  (37) 

AMINF,  ALFAD,  REREAL, 

TIME 

2546 

WRITE  (37) 

(  (  01 (I ,K) ,  1-1 ,IKAX 

),  K-l, KMAX  ) 

2547 

WRITE  (37) 

(  (  02(1 ,K) ,  I -1, IMAX 

),  K-l, KMAX  ) 

2S48 

WRITE  (37) 

(  (  Q3< I ,K) ,  1-1 , IMAX 

)  ,  K*  1  ,  KMAX  ) 

2549 

WRITE  (37) 

(  (  04(1 ,K) ,  1*1, IMAX 

),  K-l , KMAX  ) 

2550 

REWIND  37 

2551 

END  IF 

2S52 

I F  ( 

ITN  .EQ  8 *NOUT  )  THEN 

2553 

REWIND  38 

2554 

WRITE  (38) 

IMAX  ,  KMAX 

2555 

WRITE  (38) 

(  (  X ( 1 ,K) ,  1*1. IMAX 

,  K* 1 , KMAX  ) 

2556 

WRITE  (38) 

<  (  Z(I ,r ) ,  1-1, XMAX 

),  K-l, KMAX  > 

2557 

WRITE  (38, 

IMAX  ,  KMAX 

2558 

WRITE  (38) 

AMINF,  ALFAD,  REREAL, 

TIME 

2559 

WRITE  (38) 

<  (  01 ( I  , K ) ,  J*1  , I  MAX 

).  K-l.  KMAX  ) 

2560 

WRITE  (38) 

(  (  02  (  I  ,  F  )  .  1  *  1  ,  I  MAX 

),  F- i . KMAX  ) 

2561 

WRITE  (38) 

(  (  03 ( I , F ) ,  1-1 .IMAX 

),  F  *  1 , KMAX  ) 

2562 

WRITE  (38) 

(  (  Q4(l ,F),  1*1,2 MAX 

),  X-l.FMAX  ) 

2563 

REWIND  38 

2564 

END  IF 

2565 

IF  < 

ITN  -EQ  9  *NOL'T  )  THEN 

2566 

REWIND  3V 

2567 

WRITE  ()9, 

I  MAX  ,  KMAX 

2566 

V,  ’Tt  ( JO  ; 

(  (  X ( I , F ) .  I "I , IMAX 

.  F - I . FMAX  ) 

2569 

Nk.TE  (39) 

(  (  Z(  I  ,  K  )  ,  1*1  ,  IMAX 

).  F  -  1  ,  KMAX  ) 

2570 

WRITE  (39 

IMAX  .  KMAX 

2571 

W?  I  TF  (39 

AMINF.  ALFAD.  RFPf  A I. . 

TI*I 

SOURCE  PROCRAM 

nse2d.f 

msn 


18  -r. ? 


DATE 


1/20/90 


TIME  1 1 :46.-04  am 


PACE* 


36 


UNE* 


SOURCE  TEXT 


2S72  i 
^573  ! 
“2574 
f“2S7S  i 
“2576 
,  2577  j 
“2578  ' 
“2579  ■ 
“2580 
“2581 

2582 

2583 

2584 

2585 

2586  ; 

2587 

2588  i 
"2589  I 

,2590  | 
“259  j 
“2  592 
125931 
“2594  i 

2595 

2596  j 

2597  I 
i  2598  i 
1^2599  j 

2600  1 
2601  i 
'2602  1 

2603 

2604 

2605  1 

2606 

2607  ' 

2608 
2609 
“2610 
261  I 

“2612 

2613 

“2614 

2615 

2616 

2617 

2618 

2619 

2620 
2621 
2622 

2623  i 

2624 

2625 

2626 
“26 27 

2628 

2629 


WRITE  (33) 
WRITE  <J»> 


(  (  01(1, R),  1-1, IMA*  ) 


*  %  _  _ _  .  .  I« 

RRiix  ( Ji  i  (  (  Q2(I  >K)  i  ),  1“ 

WRITE  (35)  {  (  Q3 (I ,K) ,  I-l.IMAX  ),  E- 

WRITE  (35)  (  (  Q4 ( I , K )  ,  I-l.IMAX  ),  X- 


EKZ>  IF 


EMAX  ) 
EMAX  ) 
EMAX  ) 
EMAI  ) 


IF(  ITS  .EQ  10*NOUT  )  THEN 
REWIND  40 

WRITE  (40)  I MAX  ,  EMAX 
WRITE  (40)  (  (  X(1,K)»  I-1,IMAX  ),  K- 

WRITE  (40)  (  (  2(1, E).  I-l.IMAX  ),  K1 

WRITE  (40)  I MAX  ,  EMAX 

WRITE  (40)  AMINF,  ALFAD,  REREAL,  TIKE 
WRITE  (40)  (  (  Q1(I,I),  1*1, IMAX  ),  X 

WRITE  (40)  (  (  Q2(I,E>,  I-l.IHAX  >,  K- 

WRITE  (40)  (  <  Q3  ( I , K) ,  1*1, IMAX  ),  I!' 

WRITE  (40)  (  (  04(1,1),  1*1, IMAX  ),  K' 

END  IF 

ir<  ITN  .E0.  11 *NOUT  )  THEN 
REWIND  41 

WRITE  (41)  IMAX  ,  EMAX 

WRITE  (41)  (  (  2(1, X),  1*1 , IMAX  ),  K 

WRITE  (41)  (  (  1(1,1),  I-l.IMAX  ),  E 

WRITE  (41)  IMAX  ,  IMAX 

WRITE  (41)  AMINF,  ALFAD,  REREAL,  TIKE 

WRITE  (41)  <  (  01(1, K),  1-1, IMAX  ),  K 

WRITE  (41)  (  (  Q2 ( X , K ) ,  1*1 , IMAX  ),  K 

WRITE  (41)  (  (  03(1,2),  1*1, IMAX  ),  X 

WRITE  (41)  (  (  Q4(l, K),  1*1 , IMAX  ),  X 


EMAX 

EMAX 


END  IF 


If(  ITN  .EQ.  12«NOOT  )  THEN 
REWIND  42 

WRITE  (42)  IMAX  ,  EMAX 
WRITE  (42)  (  (  X( I , K ) ,  1-1, IMAX  ) 

WRITE  (42)  (  (  Z( I ,E ) ,  1*1, IMAX  ) 

WRITE  (42)  IMAX  ,  EMAX 
WRITE  (42)  AMINF,  ALFAD 

WRITE  (42)  <  (  01(1, K), 

WRITE  (42)  (  <  02 ( I , E ) , 

WRITE  (42)  (  (  03(1, K), 

WRITE  (42)  (  (  Q4 (I,K) , 

END  IF 

I F(  ITN  .EQ.  13 *NOUT  )  THEN 
REWIND  43 

WRITE  (43)  IMAX  ,  RMAX 

WRITE  (43)  (  (  2(1.2),  1*1, IMAX  ),  X 

WRITE  (43)  (  (  2(1, K),  1=1, IMAX  ),  X 

WRITE  (43)  IMAX  .  RMAX 

WRITE  (43)  AMINF,  ALFAD,  REREJL,  TIKE 


REREAL,  TIKE 
1*1, IMAX  ).  X 
1*1, IMAX  ),  X 
1*1, IMAX  ),  X 
1*1, IMAX  ),  X 


WRITE  (43)  (  (  Ql(I, X),  3*1, IMAX  ) 
WRITE  (43)  (  (  02(1, X),  1*1 , IMAX  >, 
WRITE  (43)  (  (  03(1, K),  1*1, IMAX  ), 
WRITE  (43)  (  (  04(1, E),  1*1, IMAX  ), 
END  IF 

IF(  ITN  .EQ.  14 *N0L*T  )  THEN- 
REWIND  4  4 

WRITE  (44)  IMAX  .  EMAX 

WRITE  (44)  (  (  X( I , K ) ,  1*1, IMAX  ), 

WRITE  (44)  (  (  Z(I.E),  1=1, IMAX  ), 


EMAX 

EMAX 

EMAX 

EMAX 


EMAX 

EMAX 


EMAX 

EMAX 

EMAX 

EMAX 


EMAX 

EMAX 


XMAX 
KMAX 
,  XMAX 
EMAX 


XMAX 

EMAX 


,  EMAX 
,  EMAX 
,  EMAX 
,  EMAX 


KMAX  ) 
EMAX  ) 


2630  WRITE  (44)  IMAX  ,  KMAX 

263)  WRITE  (44)  AMINF.  ALFAD,  REREAL,  TIME 

2632  WRITE  (44)  (  (  01(1, X),  1*1, IMAX  ),  K-l.RKAX  ) 

2633  WRITE  (44)  (  (  02(1, X),  1*1, IMAX  ),  X=1.KKAX  ) 

2634  WRITE  (44!  {  (  Q3(1,K),  1*1, IMAX  ),  X=1,XKAX  ) 

2635  WRITE  (44)  (  (  Q4(I, E),  1*1, IMAX  ).  K=1,KMAX  ) 

2636  END  IT 

2637  RETURN 

2638  END 


i 


ogra 


SOURCE  PROGRAM 

airfgr.f 


date  1/20/90  '’ace* 
"time  1 1:45:57  am  ^ 


SOURCE  TEXT 


PROGRAM  AIRFGRID 
PARAMETER  ( IX-300 , KX-100 ) 

C0MH0N/GR1 Dl/X { IX , EX) ,Z(IX,KX) 

COmON/OGB  I D/DT ,  I  MAI ,  KMAX , I TEL , I  TED 
LOGICAL  VISCOUS 
READ( 5,10) 

R£AD( 5,100 )  I MAX , KMAX 
R£AD( 5,10 } 

READ( 5,100)  ITEL,  ITEU 
READ( 5,10) 

READ* 5,200)  VISCOUS 
READ( 5,10) 

READ( S,300)  DM IN 
*EAD<5,10} 

READ (5,*)  AORAT  ,  AOEXP  ,  sdispl 
ILE  -  (  ITEU  -*  ITEL  )  /  2 
I  UP  -  (  ITEU  -  ITEL  )  /  2 

WRITE  (6, 1000)  IMAX ,  KMAX  ,  ITEL ,  ITEU , ILE , IUP ,  DMIN ,  AORAT 
1000  FORMAT C Ib«x  -  ‘  ,  110 , lOx , ‘ taax  -  ',110, 5x,/, 

1  'Trailing  edge  lower  -’,110,/, 

2  'Trailing  edge  upper  **,110,/, 

2  'Leading  edge  *',110,/, 

2  'lupper  *  Hover  -',110,/, 

3  ‘Distanse  of  the  first  point  -*,120.10,/, 

4  'Streching  ratio  -*,120.10,/) 

CALL  AIRFOL(  AORAT  ,  AOEXP  ,  sdispl  ) 

IF(  VISCOUS  )  CALL  CLCSTR(  DMIN  ) 

C 

WRITE( 6 , 1100 ) 

1100  FORMAT(//, 'GRID  BOUNDARIES • ,/) 

RTEOUT  -  ABS (  X(l,l)  -  X(ITEL,1)  ) 

FLEIN  -  A8S<  X< ILE, KMAX)  -  X(ILE,1)  ) 

I  UP  *  1L£  •*  <  ITEU  -  ILE  )  /  2 

RUP  -  ABS (  Z( IUP , KMAX )  -  Z(IUP,1)  ) 

WRITE (6, 1200)  RTEOUT  ,  FLEIN  ,  RUP 
120 0  FORMAT( 5X . ■ Didtaose  between  trailing  edge  ans  outflow  -* ,120.10,/ 

1  ,5x, ■ Distanse  "  leading  edge  and  inflow  -*,120.10,/ 

2  , 5x , ‘ Distanse  of  the  body  from  the  upper  boudary* ’ , f20 . 10,/) 
REWIND  21 

WRITE  (21)  I MAX , KMAX 

WRITE  (21)  ((  X(I,K),  I - 1 , IMAX ) ,  K-l.KMAX  ), 

1  <<  Z(1,T),  I-l.JMAX),  R-l ,KMAX  ) 

STOP 

10  FORMAT ( IX ' 

100  FORMAT (215) 

200  FORMAT ( 3 L? ) 

300  FORMAT ( 4 F 10 . 0 ) 

END 


1/20/90  fACt# 


TIME  11:45:57  am 


SOURCE  TEXT 


SUBROUTINE  METRIC 


SUBROUTINE  METRIC 


PARAMETER  < IX-300 , XX-100 ) 

COMMON/ FIX/OMEGA , HDOT 

COKMON/DGR I D/DT , I MAX , XMAX , I TEL , I TEU 

C0MH0N/GRID1/X<IX,XX)  ,2(11,0) 

COMHON/GR I D/Y ACOB ( I X , KX ) 

COHMON/MTRIX/XIX ( IX , EX ) , XXI < IX , EX ) , ZETAX ( IX , IX ) , ZETAZ ( IX , XX ) , 

1X1T <  IX , XX ) , ZETAT ( IX , XX ) 

C 

C***  TIE  BUB  ROUTINE  METRIC  COMPUTES  THE  METRICS  IN  ALL  THE  TWO  DIRECTIONS  AMO 
C  TIE  UNSTEADY  COEFFICIENTS  ETAT  ETC. 


DO  1000  X  -  1  ,  XMAX 
DO  1000  I  -  1  ,  IMAX 
XTAU  -  OMEGA  *  Z<I,X) 

YTAU  -  OKLGA  *  (-X(I,X))-  HDOT 

PRESEMT  SET  UP  IS  FOR  FLO»  PAST  AM  AIRFOIL. 


IP(I .CO. 1. OR. I. EQ. IMAX)  GO  TO  10 
XXI  -  .5  -  (X{I+l,f >-X(I-l,X) ) 
ZXI  •  .5  •  <Z<I+1,K)-Z<Z-2,X)) 

GO  TO  15 

10  ir<I .EQ.IKAX)  GO  TO  16 

XXI  -  1.0  •  (  X( 2 , X )  -  1(1,1) ) 

ZX1  -  1.0  «  ( Z{ 2 , X)  -  Z< X , K) ) 


ZX1  -  10  «  (  Z{ 2  ,  X )  -  Z(1,X)) 

GO  TO  15 

16  XXI  -  1.0  «  ( X < IMAX , X J  -  X(IMAX-1 ,X) ) 

ZXI  -  1.0  *  ( Z( IMAX , X )  -  Z( IMAX- 1 , X ) ) 

15  CONTINUE 

IF( X . EQ. 1. OR . X.EQ.XMAX)  GO  TO  1? 

XZET  -  .5  •<X(1,X+1)-X(I,X-1)) 

Z2ET  -  .5  •{!(! 

GO  TO  20 

17  IF(K.EQ.XKAX)  GO  TO  18 

XZET  -  2.  *  X( I , 2 )-l . 5  *  X(I,1)  -  .5  *  1(1,3) 

ZZET  »  2.  *  Z( I , 2 )  -  1.5  •  Z(I,1>  -  .5  «  Z<I,3 

GO  TO  20 

18  XZET  -  1.5  *  X ( I ,XKAX)'2.«  X( I , XMAX-1 )♦ . 5*X ( 1 , ' 
ZZET  ■  1.5  •  Z<I ,XKAX)-2.«  Z ( I , XKAX-1 )+ . 5*Z ( 1 , i 
20  CONTINUE 

YACOBI  “  XXI  •  ZZET  -  XZET  *  ZXI 
YACOB( I f  X  J  *  1.  /  YACOBI 
X1X(I,X)  -  ZZET  •  YACOB(I.X) 

XIZ(I.X)  -  -XZET  *  YACOB( I / X) 

XTAU  -  OMEGA  «  Z<I,X) 

YTAU  »  -  OMEGA  •  X(I,X)-  HDOT 

XIT( 1 , X )  -  -  XIX ( I , X )  *  XTAU  -  XIZ(I,X)  *  YTAU 
ZETAX ( I »  X )  *  -ZXI  *  Y ACOB ( I , X ) 

ZETAZ ( I  ,  X  )  -  XXI  «  YACOB<I,X) 

ZETAT ( I , X )  -  -  ZETAX ( J , X)  *  XTAU  -  ZETAZ ( I, X) 
1000  CONTINUE 
RETURN- 
END 


X(I ,XMAX-1)+.5«X(I ,KMAX-2) 
Z(I ,XMAX-1)4 . 5*Z( 1 , XMAX-2 ) 


SOURCE  TEXT 


107 

108  c«* 
T69  c* 

no  c* 

lit  c* 
Tli  c*« 
113 
IT  4 

115  t  c 
”1  16  !  c 
117  c 
TI  8  j  c 

119  c 

120 


SUBROUTINE  MRAPdl  ,JJ  ,ISING,VSING,XP, TP, SO,  AO.YO) 


SUBROUTINE  WRAP 


.124  !  1 

i  s'] 

. "726 

"127  1 
T28  i 
T29 
T30 

"-nil 

18  j 

135  { 

136 
T37 

38  i 

139  ! 

140  ] 

141  ! 

142  i  2 

1 43  ]  c 

. . i'44  i  c 

. 145  j 

146  ;  1000 

147  2000 

148 


PARAMETER  ( IX-300 , KX-100 ) 

DIMENSION  SO(IX,4),YO(IX,4) ,A0(IX,4) ,XP(1) ,YP<1) 

BUS  SUBROUTINE  UNWRAPS  T&E  AIRFOIL  AND  STORES  THE  UNWRAPPED 
SURFACE  IK  ARRAYS  AO  AND  SO.  IT  ALSO  DETERMINES  THE  STRETCHING 
IN  THE  ETA  DIRECTION. 

IMID  -  (II  +  1)  /  2 
DY  -  -8  /  (JJ  -  2) 

DO  1  J  -  2  ,  JJ 

Y  -  FLOAT (J-2 )  *  DY 

Y0( J , 1 )  -  1.25  *  Y  /  (1.  -  Y  *  Y) 

Y0( 1  ,  1)  -  -  Y0( 3 , 1 ) 

PI  -  4 .  •  ATAN  (  1 . ) 

ANGL  *  PI  +  PI 
U  -  XP(1)  -  XSING 

V  -  TP(1)  ~  liiiiG 


PI  -  4 .  •  ATAN  (  1 . ) 

ANGL  •  PI  t  PI 
U  «  XP(1)  -  XSING 

V  -  YP(1J  ~  iSiiJG 

U  -  l. 

V  -  0 

IIM1  -  IX  -  1 
DO  2  I  -  1  ,  II 
XII  -  XP(I)  -  XSING 
Yll  -  YP( I )  -  YSING 

ANGL  -  ANGL  «■  ATAN2 ( (U*Y11-V«X11 ) , (U*X11+V*Y11 ) ) 

R  -  SORT ( 111 • *2  Yll  *  *2  ) 

U  «  Xll 

V  -Yll 

R  -  SQRT(R) 

A0 ( I , 1 )  -  R  *  COS ( . 5  *  ANGL) 

S0(1,1)  -  R  •  SIN ( .5  *  ANGL) 

WRITE  (4,1000) 

WRITE  (6 , 2000}  (I, AO(I,2 ),S0(I,1 ),I  •  1,  II) 

RETURN 

FORMAT (IX, ’UNWRAPPED  COORDINATES  IN  THE  TRANSFORMED  PLANE’) 
FORMAT ( 1 5  ,  2F20.8) 

END 


SOURCE  PROCRAM 

airfgr.f 


date  1/20/90  PACE  * 

1 1X5:57  am  5 


LINE# 


2IIJC 
12 1 21 
213  jc 


SOURCE  TEXT 


SOURCE  PROCRAM 

I23HEE3S3I 

airfgr.f 

™«  11:45:57  am 

SOURCE  TEXT 


SUBROUTINE  SINO(N2,N,X,Z,SLE,YLE,TEA,TrS,XSING,rSINa,CBD| 


i*4 
85  C 
286  C 
1387  c 
1288  c 
289 
'  290 

42 

4 

m: 

16 

__29? 

>6 
►9 


SUBROUTINE  SING 


PARAMETER  < IX-300 , KX-100 ) 


tilt  SUBROUTINE  COMPUTES  SINGULAR  POINT  LOCATIONS. 

DIMENSION  X( 2 )  ,  1(2) 

NU  -  N2 
N1  •  N2  ■*  1 
M3  -  N2  -  1 
XI  -  X  ( K 1 } 

21  -  2<N1) 

12  -  X ( K  2 ) 

12  -  2(N2  ) 

X3  -  X(N3) 

S3  -  Z( N 3  ) 

D1  -  X2  •«  2  -  XI  **  2 

D2  -  22  «*  2  -  Z1  **  2 

D3  -  12  -  XI 

D4  -22-21 
DS  -  X3  **  2  -  XI  •*  2 

Dt  «  13  *«  2  -  21  *•  2 

D7  -  X3  -  XI 

DS  -  23  -  21 

B  -  (D7  *  (  D1  +  D2)  -  D3MD5+D6)  )/<  2 .  *  <D7  *D4-D3*DB  )  ) 
I P{ ABS ( 03 ) . LT . ABS ( D7 ) )  GO  TO  10 
A  ■  (D1  ♦  D2  “  2.  •  B  •  D4 )  /  ( 2 .  •  D3) 

GO  TO  20 

10  A  -  (D5  +  D6  -  2.  •  B  *  D8 )  /  (  2 .  •  D7J 
20  CONTINUE 

R  -  SQRT(  <X2-A)**  2  <Z2-B)«*2) 

XLE  -  X  ( NU ) 

YLE  -  Z{  NU  ) 

CHD  -  X(l)  -  X ( NU ) 

A2  -  (Z( 2 )~Z(1 ) )  / (X(2)  -  X(1J) 

A1  -  (Z(N)-Z<N-1))/(X(N)-X(N-1)) 

TES  -  .5  *  (A1  ■*  A2) 

TEA  -  A2  -  Al 

TEA  -  TEA  «  57.29578 

XSING  -  (A+XLE)  /2 . 

YSING  -  (B-pYLE)  /  2. 

RETURN 

END 


SOURCE  PROGRAM 

date  1/20/90 

PACE  # 

airfgr.f 

TtME  1 1:45:57  am 

7 

SOURCE  TEXT 


SUBROUTINE  AJRFOL(  AORAT  ,  «0cxp  .sdlspl) 

C*  SUBROUTINE  AIRPOL  * 

PARAMETER  < IX-300 , XX-100 ) 

COWON/GRID1/X  ( IX ,  KX  J ,  Z(  IX ,  rX ) 

COWION/DCR I D/DT ,  IMAX ,  EMAX ,  ITEL ,  ITEU 
C04G40N/Y  SYM/ISYM 

DIMENSION  S0( IX  ,  4  )  , A0( IX , 4  )  ,  Y0(  I X  ,  4  ) , XP( IX } , YP( IX  J , 

1E(IX) , F( 3  X ) , BO ( IX ) 

DIMENSION  XL(IX) ,  XU<IX),  YL<IX),  YU(IX), 

2  XX(IX),  YY(IX) 

C 

C  DATA  (BO (I), 1*1, 32 J/l . ,1*0414,1 .0X36 ,1.1 270,1. 1715, 1.2175, 1.2651 

c  11.3145,1.3459,1.4199,1.4755,1.5349,1 .5973,1.  4436,1.  7342, 1.B099, 

c  21,  «fH,i  9799, 2. 0764, .2. 1629, 2. 3012, 2. 4341, 2. 5653,2.7597,2.9646, 

C  33 .2104, 3 . 5141 , 3 . 9019,4 .4219,5.1687 ,6.3632 ,6 . 6609/ 

C 

READ{5,*)  ISYM,  IBMAX 
IP<  ISYM  .  EQ.  1  )  THEN 
DO  101  1-1 , IBMAX 

101  READ ( 5 , * )  XU(I),  YU ( 1 ) 

ELSE 

DO  102  1-1, IBMAX 

READ( 5 , * )  XU(I),  YU ( I ) ,  YL(I) 

102  CONTINUE 
ENDIF 

IT(  ISYM  EQ.  1  )  TBEN 
DO  103  1*1, IBMAX 

103  YL( I )  -  -  YU(I) 

ENDIF 

C 

DO  1000  1-1, IBMAX 
IU  -  I  +  IBMAX 
XX(IU)  •  X t ( 1 ) 

YY(IU)  -  Yl'(I) 

IL  •  IBMAX  -  I  -»  1 
XX(I)  -  XU(  I L) 

YY(1)  *  YL(IL) 

1000  CONTINUE 
C 

IBMAX2  -  2*1 BMAX 
DO  1010  1*1,1 BMAX  2 

XP(I )  -  XX ( I  ) 

YP  ( I )  -  YY  ( I ) 

1010  CONTINUE 

FNV  -  IBMAX 

FNL  -  IBMAX 

C 

C  THIS  SUBROUTINE  GENERATES  SHEAR  PARABOLIC  C-GRID 

C  THE  FOLLOWING  SUBROUTINES  ARE  RELATED  TO  THE  GRID  GENERATION 

C  WRAP  SING 

C  TABIST  CLUSTR 

C  TAINT  STRTC 


AO (1,1)  • 

do  e  i  «  2 


<  A0< 1-1,1)  *  AORAT  ) * • AOEXP 
*  AT  AN  ( 1  .  ) 


PI  »  4  *  AT 

NU  -  FNV 
NL  •  rNL 
N  •  NU  *  NL 
IT  *  ITEL 
IE  -  2  TEV 
ILE-  <  ITEL  * 
II  -  IMAX 
JJ  -  XMAX-*  1 
IIP1  -  II  ■*  1 
1 1 Ml  •  II  -  1 
IIJJ  -  II  •  JJ 


SCALE  -  1./  (  XP(1)  -  XP {NL )  ) 

DO  3333  1-1 ,N 
XP(I)  -  XP(I)  •  SCALE 
YP  ( I )  -  YP ( I )  *  SCALE 
3333  CONTINUE 

CALL  SING  (NU  ,N  ,  XP ,  YP  ,XLE,  ZU  ,  TEA ,  TES  ,  XSING ,  YSZNG  .  CHD ) 
CALL  TABI NT ( XP , YP , XS ING , YSI NG , N ,  sd  1 spl ) 


LINE# 


SOURCE  TEXT 


„4B5 

SUBROUTINE  STRTCB<SUMDX,DX1,F1,N1,R) 

-w 

c* 

* 

488 

c* 

SU BROUTINE  STRTCS 

489 

c* 

490 

C*< 

491 

PARAMETER  (IX»300,K-100) 

492 

c 

493 

c 

IBIS  SUBROUTINE  DETERMINES  A  GEOMETRIC 

494 

c 

PROGRESSION  OT  GRID  SPACING  BETWEEN  I  AND  HI  SUB  TEAT 

495 

c 

SOME  DR)  EQUALS  SUKDX .  THE  RATIO  BETWEEN  SUCCESSIVE 

496 

c 

STAGINGS  IS  R. 

'497 

N  -  Ml  -  1 

498 

ft  -  1.5 

499 

El  *  l.E-04 

400 

E2  -  l.E-04 

501 

DO  10  L  -  1 ,  50 

302 

F-  <R-1)  •  SUMDX  -  DX1M*"N-1) 

M3 

rp  -  SUMDX  -  DXI  •  FLOAT(N)  •  R  ••  (N-l) 

304 

ft  I  TER  -  R  -  F/  FP 

505 

c 

IT(l.E-02 . IT. RJ TER. AND. RITER.LT. 1*  J  RITER  -  1. 

506 

c 

TT(1. .IT. R1TRR.AND.RITER.LT. 100. )  SITES-. 01 

S67 

ir( ABS(R-RITER) . LT .  R-El)  GO  TO  1 

508 

R  -  RITER 

509 

10  CONTINUE 

SIO 

R  -  1.0001 

5  M 

c 

DXI  -  DZTOT/ FLOAT (N1 - 1 J 

517 

RETURN 

513 

1  R-  RITER 

514 

RETURN 

515 

END 

1 


n  n  o  n  a 


SOURCE  PROCRAM 


DATE 


1/20/90 


vm  11*5:57  am 


PACE# 


10 


SOURCE  TEXT 


PARAMETER  <IX-300,rX-100> 

ROTATE  GRID  IK  THE  CMCEWSE  DIRECTION  BT  AN  AMOUNT  DALTA 
COtnON/GRIDl/X(IX,n),t(IX,KZ) 

CA  •  COS (DAL FA) 

SA  -  -  SIN (DALFA) 


DO  20  K  -  1 
DO  20  I  -  1 
XI  -  X(X,K) 
IX  -  Z(I,X) 
X( I , K)  -  XI 
20  Z(I,K)  *  Z1 
RETURN 
END 


KHAX 
I  MAX 


CA  -  ZX 
CA  XI 


SA 

SA 
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